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The  major  goal  of  this  research  is  to  develop  a 
systematic  methodology  for  analyzing  the  dynamical 
properties  of  electroencephalogram  (EEG)  and  simulate  its 
dynamics  with  artificial  neural  network  models.  The  study 
assumes  that  the  EEG  is  produced  by  a nonlinear  dynamical 
system  and  promotes  the  understanding  of  complex  brain 
activity  through  network  simulation. 

This  research  starts  with  the  estimation  of  correlation 
dimension  and  the  largest  Lyapunov  exponent  of  EEG  signals. 
The  experimental  procedures  and  selection  of  implementing 
parameters  are  discussed  in  detail.  Our  results  show  that 
the  EEG  dynamics  have  a low  correlation  dimension,  which  is 
an  indication  for  an  EEG  chaotic  attractor.  But  more 
importantly,  we  developed  the  methodology  to  produce  more 
reliable  estimates  of  dynamical  characteristics  from  time 
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series.  The  effect  of  signal  condition  (filtering)  and 
time-varying  system  parameters  were  also  investigated.  The 
nonlinear  dynamics  of  EEG  segments  belonging  to  different 
brain  states  of  two  subjects  were  analyzed  to  show  the 
capability  of  our  approach  in  categorizing  the  EEG  signals. 
This  methodology  was  also  applied  to  explore  the  potential 
for  modelling  the  nonlinear  dynamical  behavior  of  the  EEG 
with  a class  of  recurrent  artificial  neural  networks.  The 
network  simulation  showed  that  a periodic  weight  variation 
in  this  neural  network  produced  an  output  that  is  visually 
similar  to  the  EEG  during  a grand  mal  seizure  or  during  slow 
wave  sleep.  Moreover  the  correlation  dimension  and  largest 
Lyapunov  exponent  estimated  for  the  output  signal  produced 
values  in  the  EEG  range. 

This  modelling  approach  can  be  used  to  help 
neurophysiologists  understand  the  effects  of  patterns  of 
activation  on  neural  assemblies  and  their  relationships  with 
measured  EEG  signals. 
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CHAPTER  1 
BACKGROUND 

Introduction 

For  decades,  the  electrical  activity  of  the  human  brain 
(electroencephalogram  - EEG)  has  been  extensively  studied  in 
the  hope  of  understanding  the  complex  nature  of  brain 
function  and  helping  clinicians  diagnose  and  treat  brain 
disfunctions.  To  deal  with  an  enormous  amount  of  recorded 
EEG,  data  guantif ication  has  been  the  primary  step  of  EEG 
analysis.  A number  of  methods,  generally  classified  as  time 
and  frequency  domain  methods,  have  been  developed  to 
quantify  the  EEG.  For  instance,  the  heuristic  method  (1,2) 
(time  domain  analysis)  serves  as  a valuable  clinical  tool 
for  the  diagnostic  of  brain  diseases  due  to  its  good 
performance  in  EEG  feature  detection  and  pattern 
recognition.  Frequency  domain  analysis  can  be  classified  as 
parametric  and  non-parametric  analysis  (3,4).  These 
techniques  and  models  for  EEG  signal  analysis  mainly 
contribute  to  the  clinical  assessment  of  human  brain  status. 
No  further  indication  in  exploring  the  complexity  of  human 
brain  activity  is  achieved.  The  non-parametric  methods 
basically  transform  the  signal  to  the  frequency  domain 
(periodogram)  (5),  with  the  hope  that  in  the  frequency 
domain  the  signal  characteristics  become  more  visible.  One 


1 


2 


of  the  problems  of  the  frequency  domain  methods  is  that  when 
the  time  domain  signals  are  of  complex  shape,  an  infinite 
number  of  sine  waves  is  necessary  to  actually  represent  the 
time  signals  in  the  frequency  domain.  In  the  stochastic 
model  (parametric  analysis,  they  are  mostly  linear  models) 
the  assumption  that  the  EEG  is  the  output  of  a linear  system 
excited  with  a random  input  actually  tells  us  that  the  EEG 
contains  no  information,  which  contradicts  the  clinical 
observation  and  cannot  explain  the  relationship  between  EEG 
and  the  brain  function.  Moreover,  there  is  no 
neurophysiological  evidence  to  corroborate  the  linearity 
hypothesis.  In  fact,  nonlinear  properties  have  been 
observed  from  the  little  that  is  known  about  the  brain  and 
the  neuron.  The  characteristics  of  massive  interconnection 
and  parallel  processing  of  the  brain  provide  also  evidence 
that  the  relationships  between  sub-systems  is  nonlinear. 

Both  the  nature  of  brain  function  and  the  problems  of 
decomposing  complex  time  domain  signals  with  (an  infinite 
number  of)  sinewaves  make  us  think  that  other  methods  of 
signal  analysis  and  quantification  may  be  more  appropriate 
to  study  the  EEG  as  an  observable  of  a highly  complex 
nonlinear  system.  Naturally,  the  recent  advances  in  the 
field  of  nonlinear  dynamics,  both  in  mathematics  and  in 
physics,  caught  our  attention. 

A new  field  of  research  in  EEG  studies  which  models  the 
human  brain  as  a nonlinear  dynamical  system  provides  great 
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illumination  for  the  investigation  of  the  macroscopic 
properties  of  brain  activity  (6-10) . It  gives  a new  way  to 
explain  the  EEG  waveforms.  Besides,  the  modelling  of  EEG 
nonlinear  dynamics  offers  a possibility  of  understanding  the 
linkage  between  the  recorded  EEG  and  brain  function,  which 
is  inaccessible  by  the  other  methods.  For  example,  the 
capability  of  the  brain  to  generate  and  process  information 
can  be  described  as  the  result  of  differences  in  nonlinear 
system  dynamics.  Accordingly,  lower  dimensionality  of  EEG 
activity  was  obtained  (8,9)  for  the  type  of  petit  mal 
epileptic  seizure  when  compared  to  an  active  mental  state 
(arithmetic  calculation) . The  agents  producing  the  seizure 
tend  to  drive  the  brain  activity  towards  a stable  periodic 
motion  (like  a limit  cycle  attractor)  so  that  information 
processing  would  be  reduced  a lot  in  such  states.  It  was 
reported  that  the  dimensionality  of  the  reconstructed  EEG 
attractor  (explained  below)  increases  with  the  increasing 
mental  activity  (6-9) . 

The  human  brain  is  probably  the  most  complex  system  we 
know;  according  to  the  present  estimation,  there  are  about 
101Q  to  1011  elements  (neurons)  in  the  hjiman_brain  and  103 
or  104  functioning  synapses  per  neuron!  It  is  absolutely 
impossible  to  cope  with  a complex  system  with  such  an 
enormous  amount  of  elements  and  interconnections  unless  some 
global  measures  of  activity  are  developed.  Fortunately, 
scientists  are  pursuing  a concept  which,  instead  of  getting 
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the  information  of  individual  unit,  deals  with  the 
macroscopic  properties  of  a complex  system  (11) . This 
procedure  can  yield  substantial  data  reduction.  Quite 
evidently  an  enormous  compression  of  information  takes 
place,  by  condensing  1010  to  1011  variables  of  the 
microworld  into  only  several  variables  of  the  macroworld. 

One  of  the  issues  is  the  selection  of  the  relevant 
quantities  to  model  the  brain  states.  The  noticeable 
progress  in  the  theory  of  nonlinear  dynamics  during  the  past 
ten  years  (12,13)  advances  proposals  for  modelling  the  brain 
as  a nonlinear  dynamical  system  having  the  property  of 
continuous  and  spontaneous  changing  state  (14) . Recent  work 
(15-19)  further  helps  one  develop  the  ability  to  investigate 
real  dynamical  phenomena  and  dynamical  systems  by  analyzing 
the  experimental  data  (20,21).  Nonlinear  dynamics,  a new 
interdisciplinary  field,  has  been  applied  not  only  to 
problems  in  physics  but  also  to  a wide  variety  of  nonlinear 
problems  in  other  areas,  such  as  the  feedback  control  of 
electrical  circuits  (22) , the  pattern  formation  and 
recognition  (23) , information  theory  (24-26) , the  evolution 
of  chemical  reactions  (22),  etc.  From  an  engineer's 
aspect,  these  new  concepts  offer  a fresh  way  to  proceed  with 
old  data  and  give  new  explanations.  In  fact,  it  has  been 
reported  (6-10)  that  the  brain  dynamics  can  be  described  by 
a low-dimensional  chaotic  system  with  deterministic  nature. 
This  observation,  which  confirms  that  the  complex  brain 
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dynamics  involves  only  a few  degrees  of  freedom,  arouses  our 
attention  of  investigating  in  depth  the  brain  modelling 
using  dynamical  theory.  The  idea  and  scheme  are  described 
in  the  following. 

First  of  all,  from  the  one-dimensional  time  series 
(e.g.  EEG)  Takens  Embedding  Theorem  (15)  enables  us  to 
generate  the  n-dimensional  signal  (the  delayed  versions  of 
the  original  time  series)  which  will  be  used  to  construct 
the  phase  space  portrait  of  EEG,  keeping  unchanged  the 
topological  invariants  of  the  system  that  produce  the 
signal.  The  selection  of  "delay"  is  an  important  issue 
since  it  affects  the  quality  of  the  reconstructed  EEG 
trajectory  and  its  quantified  dynamical  properties.  Mutual 
information  (27)  is  introduced  to  determine  this  parameter. 
The  phase  space  trajectories,  graphically  depicting  the 
dynamics  of  brain  activity,  have  been  shown  to  converge 
towards  a subset  of  the  phase  space  (i.e.  the  attractor)  for 
several  stages  of  normal  and  pathological  EEG's.  Under  the 
assumption  of  differentiable  dynamics,  this  convergence  to  a 
subspace  indicates  the  presence  of  a chaotic  attractor,  i.e. 
deterministic  dynamics  in  the  EEG.  Some  methods  derived 
from  the  nonlinear  dynamical  theory  have  been  proposed  to 
characterize  the  chaotic  attractors.  Dimensionality 
analysis  is  the  first  and  the  most  important  step  to 
quantify  the  EEG  chaotic  attractor.  In  this  study,  the 
correlation  dimension  introduced  fay  Grassberqer  and 
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Procaccia  (17,19)  has  much  appeal  for  us  due  to  its  easy 
implementation  and  fast  convergence  (with  less  amount  of 
data) . The  estimation  of  correlation  dimension  allows  us  to 
determine  the  minimum  number  of  variables  necessary  to  model 
the  dynamics  of  the  attractor  (6,8,9,28).  In  this  study,  we 
also  emphasize  the  time-varying  characteristics  of  real 
biological  systems  and  study  its  effect  on  the  correlation 
dimension. 

Our  model  to  apply  dynamical  analysis  to  the  EEG  is 
based  on  the  evidence  that  the  signal,  although  time 
varying,  has  stable  local  time  characteristics  (such  as 
sleep  stage,  seizures,  etc.).  The  system  parameters  during 
these  stationary  segments  do  not  change.  Therefore,  the 
purpose  of  the  analysis  is  to  quantify  the  dynamical 
properties  of  each  segment  with  the  tool  of  nonlinear 
dynamics. 

In  addition  to  the  dimensional  estimation  quantifying 
the  global  characteristics  of  an  EEG  attractor,  the  largest 
Lyapunov  exponent  (16,25,29,30)  is  estimated  to  probe  the 
(averaged)  local  structure  of  the  attractor.  The  spectrum 
of  Lyapunov  exponents,  the  indicator  of  the  sensitivity  of 
chaotic  dynamics  to  the  initial  conditions,  has  been  proven 
to  be  the  most  useful  dynamical  diagnostic  for  chaotic 
systems  (16).  The  approach  of  exponent  estimation  is 
basically  based  on  the  methodology  proposed  by  Wolf  et  al. 
(16)  . 
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Identification  and  characterization  of  EEG's  chaotic 
attractors  play  an  important  role  in  understanding  the 
dynamics  of  brain  activity.  In  addition,  model  construction 
has  been  a useful  tool  for  studying  the  complex  biological 
systems.  We  believe  that  network  simulation  of  EEG  dynamics 
provides  a practical  approach  for  the  quantitative  and 
qualitative  study  of  complex  brain  dynamics.  Artificial 
neural  networks  are  selected  as  the  implementation  medium 
for  the  simulation  due  to  their  nonlinear  computational 
nature,  the  massively  parallel  processing  (31,32)  and  their 
architectural  affinity  to  real  biological  neural  networks . 
Among  those  various  networks,  the  "recurrent  neural  networks 
with  a delay  function"  proposed  by  Kleinfeld  (33)  and 
Sompolinsky  and  Kanter  (34)  arouse  our  attention  and  will  be 
applied  to  the  EEG  dynamics  simulation  since  the  networks 
are  capable  of  making  transitions  between  selected  memory 
states  instead  of  resting  in  a single  state.  This  property 
has  the  potential  to  simulate  the  chaotic  signal. 

The  characterization  of  dynamical  systems  with  known 
nonlinear  differential  equation  of  motion  is  not  easy. 
However,  the  estimation  of  dynamical  properties  from  time 
series  as  in  the  case  of  the  EEG  signal  is  even  subject  to 
more  difficulties.  In  addition,  the  attractor  of  the  system 
is  not  a priori  accessible  in  many  experimental  situations 
where  neither  the  essentially  required  variables  and  the 
evolved  length  of  the  attractor  are  known.  The  approach  of 
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quantifying  the  nonlinear  dynamics  of  the  reconstructed 
trajectory  has  been  in  widespread  use  in  many  different 
cases.  Up  to  the  present,  the  published  work  (6-10)  focused 
on  the  proof  of  existence  of  EEG  chaotic  attractors,  but  no 
intensive  research  work  was  conducted  with  respect  to  the 
implementation  details  (such  as,  the  effects  of  external 
noise,  parameters  of  implementation,  number  of  data,  etc.). 
Moreover,  considerable  discrepancies  in  the  experimental 
results  (6-10)  reported  by  different  laboratories  motivate 
our  enthusiasm  for  further  investigation.  The 
inconsistencies  might  result  from  such  factors  as 
experimental  design  and  signal's  condition.  From  this 
perspective,  it  is  the  intent  of  the  research  to  explore  in 
depth  the  estimation  of  dynamical  properties  of  EEG  signals; 
and  hopefully,  the  quantitative  results  will  help  us  settle 
some  of  the  differences  of  opinion  found  in  the  literature 
and  probe  the  modelling  of  EEG  signals  (35)  as  the  output  of 
nonlinear  dynamical  systems,  such  as  the  class  of  artificial 
recurrent  neural  network  with  time  delay. 

In  this  research,  quantification  of  EEG  and  the 
network's  dynamics  conceptually  bridges  the  gap  between 
complex  brain  activity  and  the  artificial  neural  networks. 
Concretely  speaking,  if  a certain  artificial  network  model 
possesses  an  attractor  with  the  dynamical  characteristics 
closely  similar  to  the  observed  phase  trajectory  of  EEG,  we 
believe  the  neural  network  system  can  be  used  to 
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characterize  some  of  the  dynamical  properties  of  the  system 
that  produce  the  EEG  (i.e.  the  brain  state).  The  parameters 
that  we  have  in  mind  are  presently  restricted.  The 
correlation  dimension  and  the  largest  Lyapunov  exponent 
quantifying  the  "degree  of  chaos"  are  expected  to  be  able  to 
achieve  classification  of  brain  states,  which  have  clinical 
significance.  For  instance,  the  increase  in  correlation 
dimension  or  the  decrease  in  Lyapunov  exponent  may  be 
diagnosed  as  the  increase  in  mental  activity.  In  addition, 
time-varying  characteristic  is  an  important  factor  in 
nonlinear  dynamics.  This  may  drive  a network  having  cyclic 
output  pattern  to  chaos.  It  is  the  goal  of  this  research  to 
first  develop  the  methodology  for  studying  EEG  dynamics  and 
further  investigate  the  hypotheses  described  above.  We  will 
be  using  several  EEG  segments  from  different  brain  states 
belonging  to  two  subjects  for  this  work.  This  amount  of 
data  seems  reasonable  to  address  the  algorithmic 
implementation . 

In  summary,  this  research  study  involves  the 
determination  of  the  dimensionality  of  EEG  attractors,  a 
further  characterization  of  their  chaotic  properties  through 
Lyapunov  exponents,  and  the  simulation  of  EEG  dynamics  with 
artificial  neural  networks.  The  background  concepts  for 
each  one  of  these  topics  is  presented  in  the  following. 
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Concepts  and  Definitions 
Deterministic  Dynamics  and  Chaotic  Attractor 

It  was  believed  that  systems  described  by  low  order 
differential  equations  behave  in  simple  ways.  As  long  as 
these  systems  could  be  reduced  to  a few  perfectly  understood 
and  deterministic  laws,  their  long-term  behavior  would  be 
stable  and  predictable.  Now  all  that  has  changed.  In  the 
recent  twenty  years,  physicists,  mathematicians,  biologists, 
and  astronomers  have  created  alternative  ideas.  Simple 
systems  may  give  rise  to  complex  behavior;  infinitesimal 
differences  in  input  could  quickly  become  overwhelming 
differences  in  output  (36) . This  phenomenon,  the  sensitive 
dependence  on  initial  conditions,  is  called  "chaos,”  which 
mathematically  means  the  irregular,  unpredictable  behavior 
of  nonlinear  systems  (37,38).  The  term  "deterministic"  is 
used  to  emphasize  that  the  system  dynamics  can  be  described 
by  a finite  number  of  rules. 

The  term  "dynamics"  refers  to  a system's  temporal 
evolution  which  produces  time-changing  patterns.  The 
evolution  of  a system's  state  can  be  described  and  examined 
through  the  use  of  "phase  space"  (also  called  "state  space" 
(39)),  which  gives  a way  of  turning  numbers  into  pictures. 
The  evolved  trajectory  is  constructed  in  a space  with 
coordinates  defined  by  the  state  variables  of  a system. 

Chaotic  and  random  behavior  is  an  inherent  feature  of 
many  nonlinear  systems.  By  observing  phase  space 
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trajectories,  an  astonishing  geometric  regularity  was 
observed  within  chaotic  behavior;  however,  the  random 
behavior  simply  results  in  a mess  of  trajectory  filling  the 
entire  phase  space.  This  fact  shows  that  the  evolution  laws 
of  chaotic  systems  exhibit  deterministic  characteristics. 
These  systems,  which  once  seemed  completely  intractable  from 
an  analytic  point  of  view,  possess  deterministic 
characteristics  (39-41)  and  can  now  be  understood  in  a 
geometric  or  qualitative  sense  (28,42-49). 

Given  a dissipative  dynamical  system  with  a set  of 
initial  states,  as  time  goes  on  and  the  initial  transients 
dies  away,  the  state  of  the  system  approaches  a typical  area 
of  phase  space — the  attractor.  In  other  words,  the 
dissipation  of  energy  (i.e.  the  loss  of  energy)  in  a system 
will  lead  to  the  contraction  of  the  phase  space  toward  an 
attractor.  For  dissipative  systems,  the  volume  of  the 
attractor  is  in  general  very  small  relative  to  the  size  of 
phase  space.  Figure  1-1  plots  different  categories  of 
attractors.  The  phase  trajectory  in  Figure  1-1 (a),  which 
numerically  (50,51)  simulates  a 2D  nonlinear  equations 
(52,53) 

dx/dt  = x ( 3-x-y) , (1.1) 

dy/dt  = y(x-l) , (1.2) 

shows  that  the  system's  energy  is  dissipated  so  that  the 
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Figure  1-1  Attractors. 

(a)  Fixed-point  attractor; 

(b)  Limit-cycle  attractor; 

(c)  Lorenz  chaotic  attractor; 

(d)  Rossler  chaotic  attractor. 
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outer  regions  of  high  energy  are  pulled  toward  the  center 
(the  inner  region  of  low  energy) , which  is  a fixed-point 
attractor.  Here  the  trajectory  indicates  the  motion  of  a 
system's  state  variables  in  phase  space.  Figure  1-1 (b)  is  a 
2D  projection  of  a 3D  trajectory  reconstructed  from  a 
digitized  sinusoidal  waveform  (sampling  rate  = 250  Hz) 

10  sin  ( 8irt ) + 15  sin  (207rt)  , (1.3) 

with  time  delay  = 0.048  seconds.  The  attractor,  forming  a 
closed  loop  in  the  phase  space,  is  a limit  cycle.  If  a 
system's  output  varies  periodically,  the  attractor  must  be  a 
limit  cycle  with  the  system's  state  travelling  in  a loop, 
passing  through  the  same  position  in  phase  space  repeatedly. 
Another  class  of  phase  trajectories  is  the  so  called 
"chaotic  attractor"  or  "strange  attractor."  Two  famous 
chaotic  attractors,  Lorenz  (41,54)  and  Rossler  (55),  are 
shown  in  Figure  l-l(c)  and  l-l(d).  The  remarkable  fact  of  a 
system  possessing  a chaotic  attractor  is  the 
unpredictability  of  the  future  behavior  despite  the 
deterministic  nature  of  the  governing  equations.  The  beauty 
of  the  chaotic  attractor  is  the  nonperiodicity — never 
repeating  itself  and  never  falling  into  a stable  rhythm,  yet 
the  system  state  stays  inside  a limited  sub-space  of  low 
dimensionality.  These  properties  of  infinite  long  period 
and  low  dimension  make  it  possible  to  characterize  a complex 
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real-world  phenomenon  by  a simple,  deterministic  system  of 
equations.  In  fact,  even  for  a system  like  the  human  brain 
which  is  far  too  complicated  to  be  understood,  it  was  found 
that  the  phase  trajectory  constructed  from  EEG  segments 
follows  the  kind  of  motion  of  a low-dimensional  chaotic 
attractor. 

Takens  Embedding  Theorem 

The  first  problem  of  investigating  brain  dynamics  is 
how  to  construct  chaotic  attractors  from  brain  waves.  F. 
Takens  (15)  developed  a technique  of  reconstructing  the 
phase  space  of  an  attractor  from  a time  series  of  the  system 
output,  preserving  some  of  topological  invariants  of  the 
original  system.  The  idea  is  a matter  of  embedding  the  data 
in  a phase  space  of  enough  dimensions.  Here  the  "embedding" 
represents  a smooth  map  f:  X -*•  Y,  where  X is  any  data  point 
in  a time  series  and  Y is  the  mapped  n-dimensional  point  in 
the  constructed  phase  trajectory.  The  embedding  dimension  n 
corresponds  to  the  degrees  of  freedom  of  the  constructed 
phase  space.  The  foundation  of  representing  the  oner 
dimensional  time  series  by  the  n-dimensional  phase 
trajectories  contributes  to  the  exploration  of  dynamical 
properties  of  EEGs . It  has  been  recognized  that  much  could 
be  learned  about  the  dynamical  behavior  of  a system  or 
signal  from  the  analysis  of  their  constructed  trajectories 
in  a multi-dimensional  phase  space  in  which  a single  point 
characterizes  the  entire  system  at  an  instant  of  time  (56- 
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58) . Takens  Embedding  Theorem  showed  that  a set  of 
variables  describing  the  dynamics  of  the  system  could  be 
defined  and  abstracted  from  a unique  time  series.  These 
variables  are  obtained  by  shifting  the  original  time  series 
X(t)  by  a fixed  time  delay  t (r=mT  where  m is  an  integer  and 
T is  the  interval  between  successive  samples) ; therefore,  a 
point  Xi  in  n-dimensional  phase  space  is  denoted  by  a vector 
{X(tL),  Xft-t-r),  Xft^fn-l)  r ] } . By  changing  tL,  these 

points  move  in  phase  space  along  a trajectory,  which  allows 
the  development  of  the  phase  portrait  of  a signal.  The 
embedding  theorem  ensures  that  the  attractor  is  reliably 
reconstructed  under  the  restriction  of  a sufficiently  large 
embedding  dimension  n of  the  phase  space.  The  time 
evolution  of  a system  is  depicted  as  the  trajectory  in  a 
phase  space.  The  phase  space,  which  is  the  set  of  all 
possible  states  of  a system,  stores  the  complete  knowledge 
of  the  observable  brain  dynamics  during  the  data  collection 
interval.  While  the  Taken 's  approach  seems  quite  simple  and 
straightforward,  an  adequate  time  delay  must  be  chosen.  The 
method  developed  by  Fraser  and  Swinney  (59)  will  be 
introduced  next  and  applied  to  our  study  on  the  construction 
of  the  EEG  phase  portrait. 

Time  Delay  and  Mutual  Information 

The  time  delay  r can  be  chosen  almost  arbitrarily  for  a 
noise-free  signal  of  infinite  length.  However,  the  choice 
is  critical  for  a finite  time  series.  Foreseeably,  the 
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results  of  analyzing  the  dynamical  properties  of  an  EEG 
attractor  will  be  highly  biased  by  the  selection  of  time 
delay.  The  quality  of  the  phase  portraits  actually  depends 
on  the  selection  of  t (57),  which  can  be  visualized  from  the 
demonstration  in  Figure  1-2 . The  attractors  in  Figure  1- 
2(a)  and  1-2 (b)  are  constructed  from  the  X-axis  time  series 
of  the  3-D  attractor  in  Figure  1-1 (c) , with  time  delay  of  10 
and  5 samples  respectively.  Apparently,  the  structure  is 
basically  identical  to  the  original  3-D  attractor.  However, 
the  reconstructed  attractor  in  Figure  1-2 (c)  is  highly 
distorted  because  r is  chosen  so  large  that  nr  (n:  embedding 
dimension)  is  much  larger  than  the  orbital  period  of  the 
attractor  (16),  which  causes  a highly  uncorrelated  process 
(60) . On  the  other  hand,  the  magnitude  of  t should  not  be 
chosen  too  small  since  otherwise  the  attractor  in  n-D  space 
will  be  stretched  along  the  diagonal  (Figure  1-2 (d) ) . 

The  approach  developed  by  Fraser  and  Swinney,  which 
stated  that  a minimum  in  mutual  information  (MI)  provides  a 
good  criterion  for  the  selection  of  time  delay  r (59) , will 
be  applied  to  our  study  of  the  phase  space  reconstruction 

for  the  EEG  attractor.  Mutual  information  measures  the 

^ - - — — 

general  dependence  of  two  variables,  (a  nonlinear 
measurement)  and  provides  a better  criterion  for  the  choice 
of  time  delay  than  the  autocorrelation  function  (a  linear 
measurement)  which  has  been  wide  used  before  to  determine 


delay. 
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(continued) . 
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The  concept  of  mutual  information  actually  originates 
from  communication  and  information  theory  which  concerns  the 
problem  of  noise  and  efficiency  of  channel  use.  The 
following  gives  a brief  review  of  mutual  information.  Let  S 
denote  a system  which  consists  of  a set  of  possible  symbols 
(or  data,  messages,  etc.)  slf  s2,...,  sn,  each  with  its 
probabilities  p(s1)=p1,  p(s2) =p2 , • • • > p(sn)=pn.  The  amount 
of  information  for  a symbol  s^^  is  defined  as 

ICsJ  = log2  ( 1/Pi.)  , (1-4) 

since  more  information  (or  "surprise")  will  be  obtained  if  a 
symbol  with  a low  probability  is  received.  Then  over  the 
whole  system  S the  average  amount  of  information,  which  is 
the  entropy  of  the  system  S,  is  (27) 

H (S)  = S pL  log2  (1/Pi).  (1.5) 

i 


The  entropy  function  H measures  the  amount  of  uncertainty, 
surprise,  or  information  obtained  from  the  outcome  of  some 
system.  When  the  radix  of  log  is  2,  H is  in  units  of  bits. 
To  estimate  how  X(t)  correlates  with  X(t+T) , we  suppose  that 
S and  Q have  the  elements  of  {si=X(ti)}  and  {gi=X (t^T ) } , 
respectively.  First  we  examine  the  case,  "Given  a 
particular  si;  what  will  be  the  uncertainty  of  q?" . The 
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conditional  entropy  derived  as  follows  answers  our  question, 

H (Q  | S)  = 2 p(Si)  H(Q|Si) 
i 

= 22  p(sL)  p(qj|si)  log2  (l/pfq-jlsj) 

i j 

= 22  pfs^qj)  log2  ( 1/p  (qj  | s±) ) , (1.6) 

i j 

where  p(qj|si)  is  the  conditional  probability,  p(sifqj) 
gives  the  joint  probability.  Given  an  si,  the  amount  of 
difference  between  the  information  uncertainty  before  and 
after  a measurement  of  qi  is  called  mutual  information 

I(Sj.,qi)  = log2  [l/P(Si)  ] - log2  [l/p(si|qi)  ] . 

(1.7) 


The  system  mutual  information  (i.e.,  the  mutual  information 
of  a whole  system) , which  will  be  briefly  called  mutual 
information  MI  hereunder,  can  be  derived  from  the  above 
equations  as: 

I (S , Q)  = 2 p(Si)  Ifs^Q) 
i 

= 2 p(Si)  {2  p(qj|si)  I(Si,q.)} 

=22  p(si,qj)  I ( Si , q j ) 

i j 
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= 2 S p(si,qj)  {log2  [p(si,qj)  ] 

- iog2  tP(si)] 

i j 

- log2  [p(q-j)  ]} 

= H (S)  + H(Q)  - H(S,Q) 

(1.8) 

= H (S)  - H (S  | Q) 

(1.9) 

= H (Q)  - H (Q  | S)  , 

(1.10) 

which  will  always  be  a non-negative  value.  Since  {s,q}  is 
defined  as  {X(ti) ,X(tL+T) } , Q is  a delayed  process  of  S. 
I(S,Q)  computes  the  amount  of  redundancy  of  the  second  axis 
(the  coordinate  of  Q) . It  offers  a key  to  quantifying 
spatial  coherence  and  chaotic  property.  Instead  of  only  the 
first  local  minimum  of  MI  as  proposed  in  (59) , a few 
succeeding  local  minima  will  be  also  investigated  to 
optimize  the  selection.  To  speed  up  the  calculation,  the 
data  are  normally  coarsely  quantized  by  selecting  the  "box 
size"  (59) , but  the  location  of  the  MI  minima  should  be 
unchanged. 

Correlation  Dimension 

Dimensionality  analysis  is  perhaps  the  first  and 
principal  step  for  investigating  and  characterizing  the 
dynamical  properties  of  a nonlinear  system.  Then  the  other 
measurements  of  the  system's  dynamics,  such  as  the  rate  of 
information  flow  and  the  system's  entropy,  can  be  derived 
with  the  knowledge  of  the  dimension  of  the  system.  The 
dimension  provides  a lower  bound  on  the  number  of  variables 
needed  to  model  the  system ' s dynamics.  _ 
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The  concept  of  "fractional  dimension"  was  introduced 
when  the  Euclidean  geometry  failed  to  capture  the  essence  of 
irregular  shapes.  Fractional  dimension  characterizes  the 
degree  of  roughness  or  irregularity  in  an  object.  The 
attractors  of  those  fractal  objects  are  typically  chaotic 
(45) . Several  approaches  corresponding  to  different 
definitions  have  been  derived  for  the  dimensional  evaluation 
(17,19,28,43,47,48,53,61-66).  Generally  speaking,  there  are 
two  different  definitions  of  dimension  (49,67)  which 
correspond  to  different  concepts.  The  first  concept  of 
dimension,  information  dimension  and  correlation  dimension, 

i ——i  ■ — 1 " 1 *n  ' 

relates  to  the  frequency  with  which  a certain  trajectory 
visits  different  regions  of  the  space.  It  requires  both 
measurements  of  distance  and  probability  for  the  evaluation. 
The  other  definition,  such  as  Hausdorff  dimension,  only 
requires  geometric  measurements.  The  concept  of  correlation 
dimension  measurement  is  based  on  the  computation  of  the 
probability  that  any  pair  of  points  on  the  attractor  fall 
into  the  same  cube  of  the  partitioning  (66) . Alternatively 
speaking,  it  measures  the  chance  that  two  arbitrary  points 
will  be  separated  by  a distance  < r. 

The  appeal  of  correlation  dimension  m (17,19),  a lower 
bound  to  the  Hausdorff  dimension  D and  information  dimension 
o (m  < a < D)  (49,66),  is  due  to  its  easy  implementation  and 
fast  convergence  with  less  amount  of  data  (68) . In 
particular,  it  has  been  useful  in  characterizing 
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experimental  data.  Correlation  dimension,  proposed  by 
Grassberger  and  Procaccia,  measures  the  spatial  correlation 
between  points  of  a long-time  series  on  the  constructed 
attractor;  accordingly,  it  evaluates  the  coherence  of  the 
dynamical  activity  of  the  recorded  EEG  signals.  The 
correlation  dimension  is  computed  as  the  saturation  value  of 
the  correlation  exponent  vs.  embedding  dimension  curves 
(which  will  be  abbreviated  to  "CEM  curve")  for  small  r.  The 
correlation  exponent  is  the  slope  of  the  selected  linear 
region  of  the  logarithmic  correlation  integral  curves  (which 
will  be  abbreviated  to  "CIM  curve") . Figure  1-3  displays  an 
example  of  CIM  curves  for  sleep  stage  two  EEG  (3750  samples, 
r = 9) . Consider  an  N-point  attractor  reconstructed  from  a 
long-time  series  (at  least  (N  + (n-l)r)  samples,  n:  the 
embedding  dimension),  the  correlation  integral  function  C(r) 
is  defined  as  (17,19,45,66) 

N 

C(r)  = lim  ( 1/N2 ) EE  0(r  - \XL  - X, | ) , (1.11) 

N-*»  i , j = l 

where  0(d)  is  the  Heaviside  function,  0(d)  = 1 for  d > 0 and 
0(d)  = 0 for  d < 0.  The  point  Xi  on  the  n-dimensional 
attractor  is  derived  from  the  original  time  series  X(t)  by 
calculating  {X(ti),  X(tL+r) , ...,  X[ti+(n-l) t] } (15),  where 
each  value  corresponds  to  a coordinate  of  Xi(  r and  n are 
variables  in  the  calculation. 

If  the  CEM  curve  is  saturated  (i.e.  the  correlation 


0 0 


24 


'c 


Figure  1-3  CIM  curves  for  sleep  stage  two  EEG. 
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exponent  becomes  constant)  beyond  some  relatively  small 
value  of  embedding  dimension  n,  and  the  system  dynamics  are 
differentiable,  then  the  system  that  produced  the  time 
series  possesses  an  attractor.  The  saturation  value  of 
correlation  exponent  is  the  correlation  dimension  of  the 
time  series,  while  the  embedding  dimension  n corresponds  to 
the  dimension  of  the  constructed  attractor.  Moreover,  the 
value  of  n also  sets  a lower  bound  of  the  number  of 
variables  required  to  describe  the  system's  dynamics. 
Generally  the  attractor  is  chaotic  if  the  estimated 
correlation  dimension  is  a noninteger. 

Lyapunov  Exponent 

From  the  point  of  information  theory,  sensitive 
dependence  on  initial  conditions  serves  not  to  destroy  but 
to  create  information.  For  instance,  in  a dissipative 
system  information  is  destroyed  by  the  contracting  flow  in 
phase  space.  On  the  contrary,  the  information  is  created  by 
the  expanding  flow  in  a chaotic  system  which  exhibits  the 
unpredictability — the  extreme  sensitivity  to  the  initial 
conditions  (20,25).  Therefore,  strange  attractors  may 
create  information  through  the  process  of  geometrical 
variation,  i.e.  stretching  or  folding,  expanding  or 
contracting.  Lyapunov  exponents  r (bits/sec),  the  average 
exponential  rates  of  divergence  or  convergence  of  nearby 
orbits  in  phase  space  (69) , provide  a measure  of  these 
topological  qualities  that  correspond  to  such  concepts  as 
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unpredictability,  and  characterize  the  disorderly  property 
of  a chaotic  system  (16,29,30,70).  An  n-dimensional 
dynamical  system  involves  n Lyapunov  exponents.  Considering 
the  long-term  evolution  of  an  infinitesimal  n-dimensional 
ellipsoid  of  initial  conditions,  the  ith  Lyapunov  exponent 
is  defined  by  (16,70) 

Tl  = lim  (1/t)  In  [Si(t)  / ei(0) ] , (1.12) 

t-x» 

where  ei(0)  and  ei(t)  are  the  initial  and  final  length  of 
the  ith  principal  axis.  A chaotic  attractor  must  have  at 
least  one  positive  Lyapunov  exponent.  The  magnitude  of  the 
positive  exponent  reflects  the  time  scale  on  which  system 
dynamics  become  unpredictable  (25,70).  For  example,  assume 
system  A and  B (3-D  systems)  have  the  initial  points 
specified  with  the  same  accuracy  (the  same  number  of  bits) , 
the  duration  (in  seconds)  of  predictability  tA  is  longer 
than  tB  if  the  largest  Lyapunov  exponent  rA  is  smaller  than 
rB.  The  standard  technique  for  calculating  the  complete 
Lyapunov  spectrum  has  been  well  developed  for  systems 
specified  by  dynamical  equations  (25,29,30,69).  However, 
the  procedure  cannot  be  applied  directly  to  our  EEG  signals 
due  to  the  limited  resolution  and  signal  to  noise  ratio  of 
the  data,  as  well  as  the  finite  amount  of  observation  of  the 
time  series.  The  methodology  developed  in  (16) , for  the 
first  time,  provides  the  possibility  of  estimating  the 
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largest  Lyapunov  exponent  for  the  experimental  time  series. 
In  this  study,  their  method  will  be  extended  and  applied  to 
EEG  signals.  The  algorithm  is  described  in  Appendix  as 
well.  Some  problems  and  uncertainties  encountered  in  the 
exponent  estimation  will  be  discussed  in  Chapter  3. 

Recurrent  Neural  Networks  with  Time  Delays 

Identification  and  characterization  of  EEG's  chaotic 
attractors  enable  us  to  investigate  dynamical  properties  of 
brain  activity  through  the  recorded  EEG  time  series.  Model 
construction  and  simulation  is  another  useful  tool  for 
studying  the  complex  biological  systems,  allowing  us  to 
explore  possible  biological  mechanisms  for  chaotic  activity 
(35) . Artificial  neural  networks,  originated  from  the  work 
of  some  neurologists  and  psychologists  (31) , have  been 
extensively  applied  to  several  interdisciplinary  areas.  The 
features  of  highly  interconnected  topology  and  massively 
parallel  processing  are  appealing  to  solve  complex 
computational  problems  (32)  and  to  achieve  human-like 
performance  in  the  fields  of  speech  and  image  recognition 
(71).  Most  recently  some  research  works  (33,34,72,73) 
explored  another  class  of  models.  Unlike  the  conventional 
networks  dealing  with  the  fixed-point  attractors  (stationary 
states)  (32,74,75),  these  modified  network  models  turn  the 
attractors  into  transients  and  force  the  system  dynamics  to 
move  around  in  the  space  of  all  embedded  memories.  These 
networks  have  been  designed  to  generate  cyclic  patterns  to 
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study  the  motor  neuron  activity  involved  in  many  rhythmic 
movements.  In  our  study  of  modelling  the  EEG  as  the  output 
of  a nonlinear  dynamical  system,  the  "recurrent  neural 
networks  with  a delay  function"  proposed  by  Kleinfeld  (33) 
and  Sompolinsky  and  Kanter  (34)  will  be  used  because  the 
networks  essentially  capture  some  important  aspects  observed 
in  biological  neural  network  (e.g.  parallel  processing, 
massive  interconnectivity , time  delay,  and  nonlinearity) . 

The  networks  are  modifications  of  the  Hopfield's  model 
neural  circuits  (32,76-78).  In  the  following  we  briefly 
introduce  the  structure  of  the  network. 

The  circuit  diagram  shown  in  Figure  1-4  displays  the  N- 
PE  (processing  element)  Hopfield  network  with  a delay 
function  D(t,rD).  Rectangles  at  intersections  represent 
resistive  connections  (tij  and  d^ ^ ) between  outputs  and 
inputs.  Dashed  lines  correspond  to  inverted  output  signals. 
The  blocks  symbolized  by  D(r)  consist  of  time  delays  and 
linear  buffering  amplifiers.  Assuming  that  the  network  is 
capable  of  storing  n embedded  patterns  Mv  (v  = 1,  ...,  n) , 
the  strength  of  each  synapse  is  (79) 


n 


Tij  = ( 1/N)  Z M£v  MjV  (i  f j) , 


v=l 


(1.13) 


(1/N)  Z Miv+1  MjV  (i  f j)  , 


v=l 


(1.14) 


where  m is  the  length  of  each  independent  sequence  and  T^  = 
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Figure  1-4  Circuit  diagram  of  time-delay  artificial  neural 
network . 
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Dii  = °*  With  a delay  function  D(t,rD)  = <S(t-rD),  the 
output  VjD  (t ) = Vj  (t-rD)  . Then  the  network  dynamics  can  be 
described  by  the  following  equation 

Ui(t)  + t i [dU^tJ/dt] 

N N 

= S TLj  V.(t)  + S DLj  Vj(t-TD),  (1-15) 

j=l  j=l 

where  ri  = RiCi  is  the  charging  time  of  the  ith  PE  and  Ri  is 
derived  from 

N N 

i/Rj.  = l/ri  + s (l/t^)  + s (l/d^).  (1.16) 

j=l  j=l 

The  magnitudes  of  the  weights,  Tij  and  D^,  are  related  to 
the  resistances. (tij , d^ , and  Ri)  in  the  circuit  of  Figure 
1-4  by  | Ti;j  | = Ri/tij  and  | Di:j  | = Ri/dij,  respectively. 

However,  one  problem  arises,  i.e.,  the  networks  cannot 
generate  sequences  with  diversified  patterns  [a65-66]  as 
observed  in  the  EEG.  Only  cyclic  sequences  were  produced  in 
these  networks  (33,34)  as  can  be  expected.  However,  we 
found  that  by  changing  the  value  of  the  weights  the  outputs 
of  time-delay  neural  networks  display  EEG-like  waveform 
having  varied  patterns.  The  correlation  dimension  of  the 
output  sequence  appears  to  be  in  the  EEG  range.  In  Chapter 
2,  we  will  introduce  a new  concept  "what  are  the  effects  on 
the  dynamical  properties  if  one  or  more  system's  parameters 
happen  to  change  during  the  evolution."  Can  we  increase,  by 
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some  mean,  the  "degree  of  chaos"  of  a system?  If  yes,  it  is 
expected  that  a simple  artificial  network  may  generate  the 
EEG-like  waveforms. 


CHAPTER  2 

CORRELATION  DIMENSION 
OF  EEG  TIME  SERIES 

The  technique  of  evaluating  the  correlation  dimension 
is  still  not  a well-shaped  procedure  when  applied  to  real 
signals  like  the  EEG,  though  it  possesses  the  advantages  of 
high  efficiency  and  easy  implementation  mentioned  before. 
Until  now,  most  of  the  published  work  focused  on  the 
evaluation  of  the  correlation  dimension.  It  is  known  that 
deterministic  chaotic  processes  produce  low  correlation 
dimension  estimates.  Therefore,  very  likely  low  values  of 
this  quantity  can  be  used  to  prove  the  existence  of  low- 
dimensional chaotic  attractors  (8) , but  no  intensive 
research  has  been  conducted  to  analyze  the  effect  of  the 
implementation  details  (such  as,  the  effects  of  external 
noise,  implementation  parameters,  etc.)  on  the  dimension 
estimates.  Noise  components  are  very  likely  to  cause  an 
overestimation  of  the  correlation  dimension.  Moreover,  the 
selection  of  implementation  parameters  can  be  misconducted. 
For  instance,  a "false  local  minimum"  may  appear  in  the  MI 
curves  for  some  noisy  signals,  which  usually  misleads  the 
selection  of  r.  Another  example  is  if  the  number  of  points 
is  not  enough  for  reliable  statistics,  or  if  the  embedding 
space  is  not  properly  chosen.  All  this  means  that 
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experimental  design  and  signal  conditioning  play  an 
important  role  in  obtaining  reliable  estimates.  We  are 
aware  that  this  empirical  evaluation  also  has  its 
limitations.  It  will  be  very  difficult  for  instance  to 
conclude  about  the  bias  of  the  correlation  dimension  as  an 
estimate.  But  we  feel  that  these  questions  involve  more 
fundamental  research  which  is  outside  the  realm  of 
engineering.  The  goal  of  this  chapter  is  to  determine 
adequate  parameters  for  correlation  dimension  computation. 
These  parameters  include  the  segment  duration  N (number  of 
samples)  which  gives  a convergent  dimensional  estimate  and 
the  time  delay  r (in  number  of  samples)  which  provides  a 
good  geometric  structure  of  the  constructed  attractor. 

Hence,  the  reliability  and  meaning  of  the  practical  results 
are  ensured.  Other  experimental  considerations  are 
discussed  as  well.  Then  with  a well-developed  methodology, 
the  effects  of  filtering  (noise)  will  be  investigated.  This 
study,  to  our  belief,  will  lead  to  the  establishment  of  a 
mature  implementation  process  for  the  dimension  estimation 
of  the  EEGs , and  furthermore,  a thorough  exploration  of  its 
value. 

In  the  last  section  of  this  chapter,  correlation 
dimension  of  the  Lorenz  model  having  a time-varying 
parameter  is  estimated.  The  study  is  stimulated  by  two  main 
ideas.  The  first  is  to  try  to  quantify  the  effect  of  time- 
varying  system  parameters  in  the  determination  of  the 
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correlation  dimension.  This  question  is  raised  because  real 
biological  systems  very  probably  possess  a time-varying 
mechanism,  i.e.  the  synaptic  strengths  in  biological  neural 
networks  are  time-varying.  This  is  the  principal  concept 
added  to  the  study  of  the  network  simulation  and 
construction.  In  that  section,  we  investigate  (1)  what  is 
the  effect  of  time-varying  parameter  on  system  dynamics,  and 
(2)  can  the  dimension  estimation  detect  the  "time-varying 
activity"  in  a dynamical  model. 

Experimental  Implementation  and  Considerations 

EEG  segments  of  sleep  stage  two  (250  Hz  sampling  rate, 
8-bit  ADC)  without  any  preprocessing  are  used  in  this  study. 
It  is  more  helpful  to  start  with  well-defined  data  for  our 
preliminary  work  of  methodology  development.  Therefore,  we 
selected  representative  segments  from  a five-minute  signal 
(CZ-A2)  recorded  from  one  patient.  The  parameters  for  the 
experimental  evaluation  are  specified  as  follows. 

Embedding  Dimension  n 

A group  of  curves  plotting  the  logarithm  of  the 
correlation  integral  versus  logarithm  of  distance  r (CIM 
curves) , with  embedding  dimension  n varying  from  2 to  28  (in 
steps  of  2) , is  evaluated  for  each  segment  with  duration  N 
and  a specified  r.  The  plot  of  CIM  curves  describe  the 
dependence  of  correlation  exponent  on  the  value  of  embedding 
dimension  n.  According  to  published  works  about  the 
dimension  estimation  of  sleep  stage  two  EEG  (6,8,9),  the 
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saturation  of  CEM  curves  occurs  at  a small  n (e.g.  n < 10). 
Hence  our  experiment  covers  this  range  for  obtaining  a 
reliable  estimate  of  /x.  To  speed  up  the  computation,  the 
CIM  curves  are  evaluated  for  every  other  n.  However  it  is 
not  an  adequate  procedure  for  obtaining  a more  precise 
estimate  of  correlation  dimension  and  the  number  of 
embedding  dimension.  Yet  the  determination  of 
implementation  parameters  is  not  affected. 

Computation  of  Correlation  Exponent 

The  linear  regression  method  (51)  is  applied  for 
calculating  the  correlation  exponent.  One  problem 
encountered  is  how  to  pick  a significant  linear  portion  for 
the  slope  computation.  A survey  of  previous  work  related  to 
dimension  estimates  shows  no  concrete  criterion  for  handling 
this  case.  However,  some  scattered  notes  (16,17,19,49), 
which  are  summarized  below,  help  us  understand  the  criteria 
for  the  selection  of  linear  regions  (assume  the  range  (rx, 
r2)  is  selected) . 

Distance  r < 

For  smajU  .distances  (r  < r-,)  , if  the  data  points 
deviate  from  the  power  law  (i.e.,  ln[C(r)]  and  ln(r)  have  no 
linear  relationship  when  In (r)  is  small),  these  points.,  win 
discarded . This  is  because  the  behavior  of  the  attractor 
is  not  yet  chaotic  under  this  selection  of  r and  n,  which 
makes  the  power  law  unable  to  hold  towards  small  values  of 
r.  In  addition,  for  small  r,  poor  statistics  (due  to  the 
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small  number  of  n-D  points  involved)  lead  to  a large  scatter 
of  points.  The  linear  region  may  be  extended  towards 
smaller  r by  increasing  N (i.e.  increasing  the  number  of 
points  in  the  estimation,  also  called  the  density  of  the 
attractor) , but  this  increases  tremendously  the  computation 
time. 

Distance  r > r2 

For  large  distances  (r  > r2)  , the  CIM  curves  saturate 
because  the  dynamical  extent  of  the  attractor  has  been 
totally  covered  under  r2.  Therefore,  data  points  beyond 
ln(r2)  will  not  be  considered  in  slope-calculation. 

Apparently,  the  meaningful  range  for  computing  the 
correlation  exponent  is  the  range  (rx,  r2)  . At  this  point, 
we  pick  up  this  range  visually  (a  method  for  automatically 
selecting  the  linear  region  may  be  developed  for  the  future 
study) . However  these  simple  rules  only  apply  to  the  case 
of  small  n (generally  speaking,  when  embedding  dimension 
does  not  exceed  number  of  degrees  of  freedom  of  the  time 
series) . The  CIM  curves  for  EEG  plotted  in  Figure  1-3 
exhibit  more  complicated  structure  not  discussed  in  the 
literature,  i.e.  the  CIM  curves  possess  more  than  one  linear 
region  separated  by  "knees"  when  n increases  beyond  a 
certain  value.  Then  a problem  emerges  "which  region  should 
we  use  for  the  slope  computation."  The  "knee"  phenomenon 
discussed  below  may  provide  guidance  for  the  estimate  of 
correlation  exponent. 
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Knee  phenomenon  in  CIM  curves 

If  the  EEG  attractor  is  constructed  in  a very  high 
dimensional  space,  "knees"  appear  on  the  CIM  curves,  which 
separate  two  linear  regions  of  different  slopes.  This 
brings  the  problem  of  selecting  the  appropriate  region.  In 
addition,  the  decrease  of  significant  linear  range  results 
in  an  increase  of  statistical  error  in  slope  computation. 

Two  factors  may  cause  the  appearance  of  knees  (49) . 
First,  a knee  will  appear  if  the  signal  from  a deterministic 
chaotic  system  is  contaminated  by  random  noise  of  smaller 
amplitude.  Secondly,  if  a system  involves  two  or  more 
subsystems  Ux,  U 2,  ...  with  corresponding  amplitudes  r1  << 
r2  <<  . ..,  knees  will  also  appear  separating  linear  regions. 
These  cases  will  be  discussed  below. 

The  principal  concept  of  the  first  case  is  based  on  the 
infinite  dimensionality  characteristic  of  random  white 
noise.  Unlike  the  deterministic  component  of  the  signal, 
white  noise  tends  to  fill  the  available  phase  space  (16) . 

The  slope  of  CIM  curves  for  random  noise  will  increase 
proportionally  to  the  embedding  dimension;  while  the  slope 
for  a deterministic  signal  whose  phase  trajectory  forms  a 

4 

strange  attractor  will  reach  a value  n (correlation 
dimension)  and  will  then  become  independent  of  n. 
Consequently,  with  a small  amplitude  random  noise  added  to 
the  deterministic  chaotic  signal,  the  CIM  plot  will  display 
two  approximately  linear  regions  separated  by  a knee.  This 
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can  be  realized  from  equation  (1.11)  defining  the 
correlation  integral  function  C(r) . The  function  C(r) 
calculates  the  ratio  of  number  of  pairs  (X^Xj)  that  are 
within  a prescribed  distance  r.  Distance  |Xi-Xj|  is  related 
to  the  signal's  amplitude  X(t)  by 

IXi-XJ  = sqrt{[X(ti)-X(tj)  ]2  + [X  (tj+r ) -X  ( tj+T ) 1 2 
+ ...  + [X(ti+(n-l)r)-X(tj+(n-irr) ]2}, 

(2.1) 

where  "sqrt{}''  is  the  square  root  function.  Thus  when  noise 
is  added  to  the  signal  with  a large  signal-to-noise  ratio, 
two  linear  regions  appear  in  the  CIM  curves.  For  values  of 
r above  the  knee  (the  contribution  from  the  higher  amplitude 
signal  component) , C(r)  continues  to  scale  like  rM.  For 
values  of  r below  the  knee,  C(r)  scales  like  rn  due  to  the 
random  jittering  of  the  trajectory  (19).  Accordingly,  for 
the  estimation  of  the  signal  correlation  dimension,  the 
range  of  ln(r)  beyond  a knee  must  be  chosen.  In  this  case, 
the  "knee"  phenomenon  can  also  be  used  to  establish  an  upper 
limit  for  the  number  of  degrees  of  freedom  of  the  embedded 
attractor  and  accordingly,  for  the  correlation  dimension  of 
the  signal.  In  fact,  increasing  n past  what  is  minimally 
required  for  the  embedding  has  the  effect  of  emphasizing  the 
level  of  noise  contamination  in  the  EEG  signal  dynamics 
(16).  Therefore,  the  number  of  degrees  of  freedom  must  be 
smaller  than  the  embedding  dimension  where  the  knee  in  the 
CIM  curves  develops. 
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The  second  case  treats  a system  with  complex  dynamics, 
i.e.  possessing  an  output  showing  a variety  of  waveform 
patterns.  A knee  dividing  two  linear  regions  is  reported  in 
the  CIM  curves  for  a compound  dynamical  system  formed  by  two 
noninteracting  subsystems  (49)  . The  range  r < rx  gives  the 
dimension  /u  of  the  complete  system,  while  in  the  range 
r > rx  the  CIM  curves  only  reveal  the  statistical 
information  (dimension  /x ' ) for  subsystem  2. 

For  sleep  stage  two  EEG,  it  is  reasonable  to  interpret 
the  knee  phenomenon  in  Figure  1-3  as  the  "intrusion"  of 
random  noise  component  since  there  exists  no  parallelism  in 
the  lower  portion  of  CIM  curves.  The  study  of  filtering 
effect  in  dimensional  estimation  may  provide  useful 
information  about  our  assumption.  However,  should  the 
complex  brain  dynamics  involve  several  subsystems,  how  these 
subsystems  associate  with  each  other  is  still  inaccessible 
at  this  stage. 

Time  Delay  t for  Reconstructing  the  Phase  Trajectories 

The  concept  introduced  by  Fraser  and  Swinney  (59) 
states  that  the  calculation  of  mutual  information  (MI)  for 
the  time  series  provides  an  excellent  criterion  for  choosing 
the  time  delay  t.  It  will  be  applied  in  this  study.  From 
the  MI  vs.  r curves,  we  can  appropriately  select  four 
different  values  of  r (located  at  local  minima)  to  calculate 
the  corresponding  correlation  dimension. 
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Duration  of  EEG  Signal 

Since  the  dimension  calculation  is  very  time-consuming 
(computation  time  is  proportional  to  N2) , it  is  a practical 
consideration  to  choose  the  minimum  duration  of  the  signal 
that  ensures  an  adequate  coverage  of  the  EEG  attractor 
dynamics.  Note  that  this  parameter  (duration  of  the  signal) 
is  highly  dependent  on  the  convergence  rate  of  the  EEG 
attractor.  The  correlation  dimensions  of  EEG  segments  with 
different  lengths,  including  4,  8,  15,  and  30  seconds  (i.e., 
1000,  2000,  3750,  and  7500  samples)  are  evaluated.  All  the 
segments  belong  to  the  same  state  of  brain  activity  (sleep 
stage  two) . 

Experimental  Results 

The  30-second  segment  of  sleep  stage  two  EEG  is  shown 
in  Figure  2-1.  Figure  2-2  to  2-5  display  the  dependence  of 
the  CIM  curves  for  different  number  of  samples  N = 1000, 
2000,  3750,  and  7500,  respectively.  The  variation  of  CIM 
curves  according  to  different  values  of  t is  plotted  for 
each  N.  The  straight  line  in  Figure  2-6  and  2-7  represents 
theoretical  dependence  for  a white  random  process.  In 
Figure  2-6,  the  CEM  curves  for  EEG  are  plotted  as  follows:  a 
starred  curve  (t  = 2),  a crossed  curve  (r  = 9),  a circled 
curve  (t  = 13),  and  a diamonded  curve  (r  = 16  for  N = 3750 
and  r = 17  for  the  others) . The  deviation  of  the 
correlation  exponent  (/*)  vs.  embedding  dimension  (n)  curves 
from  this  straight  line  is  significant  for  the  EEG  data, 
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Figure  2-2  CIM  curves  for  a 1000-sample  EEG  segment. 

(a)  t = 2;  (b)  t = 9;  (c)  r = 13;  (d)  r = 17. 
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Figure  2-2  (continued) . 
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Figure  2-3  CIM  curves  for  a 2000-sample  EEG  segment. 

(a)  t = 2;  (b)  t = 9;  (c)  r = 13;  (d)  t = 17. 
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Figure  2-3  (continued) . 
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Figure  2-4  CIM  curves  for  a 3750-sample  EEG  segment. 

(a)  t = 2;  (b)  r = 9;  (c)  r = 13;  (d)  r = 16. 
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Figure  2-4  (continued) 
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Figure  2-5  CIM  curves  for  a 7500-sample  EEG  segment, 
(a)  t = 2;  (b)  t = 9;  (c)  t = 13. 
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Figure  2-5  (continued) 
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CEM  curves  for  different  lengths  and  delays. 

(a)  N=1000 ; (b)  N=2000;  (c)  N=3750;  (d)  N=7500. 


Figure  2-6 
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Figure  2-6 


(continued) 
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Figure  2-7  Effect  of  EEG  duration  on  CEM  curves, 
(a)  t = 9;  (b)  t = 13. 
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which  prompts  the  interpretation  of  the  EEG  signal  as  the 
output  of  a non-linear  dynamical  system  with  a chaotic 
attractor . 

General  Characteristics 

An  evident  feature  of  the  CIM  curves  (Figure  2-2  to  2- 
5)  is  the  tendency  towards  parallelism  when  n is  increased. 
Accordingly,  pronounced  saturation  of  the  CEM  curves  (Figure 
2-6  and  2-7)  is  observed,  which  characterizes  a finite 
dimensional  deterministic  process.  The  value  of  n for  which 
the  CEM  curves  saturate  corresponds  to  the  order  of  the 
embedding  space  of  the  constructed  attractor.  According  to 
our  analysis  of  sleep  stage  two  EEGs , the  parallelism  starts 
with  n = 10.  In  the  case  of  small  n,  the  CIM  curves  clearly 
show  a large  scattering  of  points  (i.e.,  a significant 
deviation  from  the  linear  relationship)  in  the  region  of 
small  ln(r)  (the  lower  portion  of  the  curves) . However, 
this  phenomenon  shifts  to  the  middle  portion  of  the  curves 
when  the  embedding  dimension  is  increased.  As  mentioned 
before,  we  interpret  this  phenomenon  as  noise  contamination. 
The  effect  of  random  noise  component  is  emphasized  in  the 
CIM  curves  when  the  EEG  trajectory  is  constructed  in  a high 
embedding  dimension.  In  the  region  below  the  "knee",  the 
slope  continuously  increases  with  n,  which  characterizes  the 
random  process.  However,  in  the  upper  region  the  CIM  curves 
show  a chaotic  attractor  of  relatively  low  order.  The 
saturation  of  the  CIM  curves  occurs  in  the  region  above 
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ln(r)  = 4.6  ± 0.5,  i.e.,  the  process  of  covering  the  EEG 
attractor  in  extent  is  all  completed  beyond  this  region. 
Here,  saturation  means  that  ln[C(r)]  becomes  a constant, 
therefore  independent  of  ln(r) . 

Apparently,  the  minimum  of  ln(r)  required  to  confine 
any  two  points  in  a trajectory  increases  with  n; 
simultaneously,  the  length  of  linearity  decreases  as  n 
increases  (Figure  2-2  to  2-5) . This  situation  causes  much 
difficulty  in  slope-computation,  and  a large  statistical 
error  in  calculating  correlation  exponent  is  expected  due  to 
the  shrinking  of  linear  region.  The  occurrence  of 
intermediate  "knees"  in  the  CIM  curves  (at  ln(r)  =4.2  ± 

0.3)  visible  for  large  r indicates  that  the  embedding 
dimension  used  to  construct  the  EEG  attractor  is  already 
higher  than  the  system  order.  From  Figure  2-2  to  2-5,  the 
correlation  dimension  must  be  smaller  than  12,  since  a knee 
appears  when  the  embedding  dimension  is  equal  to  or  greater 
than  12 . 

Effect  of  Time  Delay 

To  conduct  this  part  of  experiment,  mutual  information 
MI  (59)  of  EEG  segments  is  calculated  first.  The  algorithm, 
which  operates  for  time  series  data  points,  is  also  helpful 
for  choosing  the  time  delay  for  the  higher-dimensional  phase 
space  reconstruction.  According  to  Fraser  and  Swinney  (59) , 
choosing  a "box  size"  has  its  advantages  and  liabilities. 

For  a given  number  of  data  points,  a large  box  size  causes 
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the  underestimation  of  MI,  whereas  a small  one  leads  to  an 
overestimation  of  MI.  Three  different  box  sizes  are  used  in 
this  study  to  investigate  the  effect  of  box  size  on  the 
calculation  of  MI.  Figure  2-8  plots  the  MI  vs.  r curves  for 
N = 7500  (solid  curve),  3750  (0.1"  dashed  curve),  and  2000 
(0.05"  dashed  curve),  with  the  centered  symbol  of  stars, 
circles,  and  diamonds  corresponding  to  the  box  size  of  2,  4, 
and  8.  The  MI  vs.  r curves  show  that  the  local  minima  occur 
at  the  same  values  of  time  delay  (such  as  t = 2,  9,  13,  and 
17)  for  different  box  sizes.  It  indicates  that  the  local 
minima  of  MI  vs.  r are  independent  of  the  box  size  selected. 
Therefore  for  our  study  the  box  size  is  not  a determining 
factor. 

For  each  N-sample  segment,  four  different  values  of  r 
(2,  9,  13,  and  17)  are  selected  to  reconstruct  the  phase 
trajectories  and  to  estimate  the  corresponding  correlation 
dimension.  The  CIM  curves  for  t = 2 (Figure  2-2 (a),  2-3 (a), 
2-4(a),  and  2-5(a))  display  no  abrupt  knees  as  n increases; 
instead,  the  curves  bend  gently  and  the  linearity  is 
blurred.  This  makes  the  effect  of  the  noise  component 
indistinguishable.  Moreover,  it  is  difficult  to  determine 
the  correlation  dimension  for  this  value  of  r . Visually 
speaking,  the  phase  portrait  for  t = 2 is  skewed  along  one 
direction  (Figure  2-9(a)),  very  much  like  the  effect  seen  in 
Figure  1-2 (d)  for  the  Lorenz  attractor.  For  r > 9,  the 
curves  show  obvious  knees  when  the  embedding  dimension 
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Figure  2-8  MI  curves  for  N=7500  (solid),  N=3750  (0.1» 
dashed),  and  N=2000  (0.05"  dashed). 
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Figure  2-9  Phase  portraits  constructed  from  a 3750-sample  EEG  segment. 

(a)  Delay  =2;  (b)  Delay  =9;  (c)  Delay  =13;  (d)  Delay  = 17. 
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n > 12.  Nearly  identical  phase  portraits  are  found  for  r = 
9,  13,  and  17  (Figure  2-9 (b),  2-9 (c)  and  2-9 (d) ) . Notice 
that,  for  all  the  four  different  values  of  N,  the 
correlation  exponent  (CEM)  curves  (Figure  2-6)  basically 
coincide  for  r > 9 (e.g.  /i  = 6.5946,  6.5760,  and  6.4966  for 
N = 3750  and  t = 9,  13,  and  16,  respectively).  Thus  the 
same  result  of  correlation  dimension  is  obtained  for  these 
three  values  of  r,  while  a gross  discrepancy  in  dimension  is 
observed  for  t = 2 (m  = 7.9553  for  N=3750) . As  a result,  we 
conclude  that  the  choice  of  t between  9 and  17  does  not 
actually  makes  any  difference  in  estimating  the  dimension  of 
sleep  stage  two  EEG. 

Effect  of  EEG  Duration 

Figure  2-7 (a)  and  2-7 (b)  display  the  CEM  curves  for 
t = 9 and  13,  respectively.  One  point  we  should  mention 
first  is  that  the  correlation  exponent  calculation  depends 
on  a very  small  linear  portion  of  the  CIM  curve  when  n is 
large,  especially  after  the  appearance  of  intermediate  knees 
(in  our  case,  n > 12),  which  may  affect  the  accuracy  of  the 
slope  estimation.  However,  in  the  CEM  curves  the  portion 
beyond  this  critical  value  of  n can  be  left  out  of 
consideration.  The  most  reliable  estimates  of  the 
correlation  dimension  are  the  values  obtained  in  the  n range 
between  the  development  of  parallelism  and  the  appearance  of 
intermediate  knee.  Inspecting  the  CEM  curves  in  Figure  2-7 
for  n < 12,  the  curve  for  N = 3750  (circles  and  solid  curve) 
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roughly  coincides  with  the  one  for  N = 7500  (diamonds  and 
solid  curve) ; while  the  curve  for  N = 2000  (crosses  and 
solid  curve)  is  close  to  the  above  two  curves  for  n < 9,  but 
inclines  towards  smaller  values  of  correlation  exponent  for 
larger  n.  Apparently,  the  curve  for  N = 1000  (stars  and 
solid  curve)  significantly  deviates  from  the  others,  which 
may  be  explained  by  the  following.  The  orbits  on  a chaotic 
attractor  are  shuffled  by  the  stretching  and  folding  process 
in  the  phase  space  (38) . The  trajectory  of  an  EEG  attractor 
is  formed  by  the  repetitive  process  of  expansion  and 
contraction.  Hence  the  dimension  analysis  of  an  EEG 
attractor  may  undergo  an  overestimation  or  underestimation 
if  the  analyzed  segment  duration  is  not  long  enough  to  cover 
the  system  dynamics  (i.e.  involve  enough  orbits  for 
averaging).  In  Figure  2-10 (a)  to  2-10 (d),  the  2D 
projections  of  a 3D  EEG  attractor  show  the  evolution  of 
phase  trajectories.  From  which  we  observe  that  some 
vigorous  oscillation  in  the  third-portion  of  the  30-second 
EEG  segment  (Figure  2-1)  results  in  the  stretching  of  the 
phase  portrait.  Therefore,  evaluating  the  first  4-second 
segment  (1000  points)  results  in  an  underestimate  of  /x  by 
about  4.5,  because  most  of  the  contracting  and  folding 
process  of  the  EEG  attractor  develops  within  this  region. 
Moreover,  the  knee  phenomenon  appears  earlier  in  the  CIM 
curves  (about  n = 8) . 
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From  Table  2-1  which  lists  the  estimated  correlation 
dimension  (by  averaging  the  correlation  exponents  obtained 
with  n = 10,  11,  and  12)  for  different  N and  r,  we  may 
conclude  that  it  is  immaterial  to  use  3750  or  7500  data 
points  for  the  estimation  of  the  correlation  dimension  of 
this  segment  of  sleep  stage  two  EEG. 

Discussion  and  Conclusion 

In  this  Chapter,  a methodology  is  developed  to  estimate 
the  correlation  dimension  of  the  EEG  embedded  attractor.  By 
using  the  cautiously  selected  parameters  (the  segment  length 
and  time  delay) , we  are  able  to  ensure  the  reliability  of 
the  results,  which  has  been  ignored  in  the  previous 
researches  (6-10) . The  conclusions  have  been  presented 
before  and  are  repeated  here  for  convenience.  We  propose 
that,  to  determine  the  N and  t for  dimensional  estimation, 
the  MI  curve  is  first  computed  for  a moderate  amount  of 
data.  The  first  local  minimum  having  the  magnitude  of  MI 
not  significantly  greater  than  the  flat  value  of  MI  is 
selected  as  the  adequate  time  delay.  For  instance,  local 
minima  at  2 and  6 samples  (Figure  2-8)  having  relatively 
large  MI  do  not  provide  reliable  estimates  and  good-quality 
attractors;  while  selecting  r between  9 and  20  samples 
actually  gives  very  close  dimensional  estimates  (Table  2-1) . 
Correlation  dimension  is  then  computed  using  the  selected 
(N, t)  set.  The  convergent  estimation  can  be  achieved  by 
repeating  the  above  procedure  with  N gradually  increased  (or 
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Table  2-1 

Correlation  Dimension  of  Sleep  Two  EEG 
with  Different  Duration  and  Time  Delay 


N 

T 

correlation 

dimension 

1000 

2 

7.4322 

+ 

0. 3887a 

1000 

9 

4.9923 

+ 

0.1431 

1000 

13 

4.6850 

+ 

0.0260 

1000 

17 

4.8003 

+ 

0.0327 

2000 

2 

7.7500 

+ 

0.2389 

2000 

9 

6.2716 

+ 

0.1926 

2000 

13 

6.1880 

+ 

0.1586 

2000 

17 

6.2138 

+ 

0.0256 

3750 

2 

7.9553 

+ 

0.6834 

3750 

9 

6.5946 

+ 

0.1409 

3750 

13 

6.5760 

+ 

0.1332 

3750 

16 

6.4966 

+ 

0.2135 

7500 

2 

7.9036 

+ 

0.5161 

7500 

9 

6.6379 

+ 

0.1277 

7500 

13 

6.6646 

+ 

0.1314 

a)  This  number  shows  the  deviation  with  respect  to  n 
(n  = 10,  11,  and  12  in  this  study) . 
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decreased) . The  delay  needs  to  be  re-evaluated  when  N 
changes . 

We  have  discussed  our  approach  of  deriving  n from 
several  estimated  correlation  exponents,  i.e.,  by  averaging 
the  correlation  exponents  in  the  range  of  embedding 
dimension  limited  by  the  value  of  n where  parallelism 
develops,  until  where  the  noise  effect  settles  in.  For  the 
sleep  two  EEG,  the  estimated  correlation  dimension  is  6.576 
± 0.133  (N  = 3750  and  r = 13),  the  average  of  correlation 
exponents  for  10  < n < 12.  While  analyzing  the  other  15- 
second  sleep  two  EEG  segments  from  the  same  subject,  the 
estimated  correlation  dimension  did  not  change  (Figure  2- 
11).  Notice  that  our  result  (n  = 10,  /x  = 6.576)  is  obtained 
from  the  noisy  signal.  Comparing  with  results  reported  by 
other  authors  (6,8),  which  showed  an  embedding  dimension  of 
6 and  the  correlation  dimension  /x  of  5.00  - 5.03,  the 
difference  may  be  caused  by  two  factors.  One  is  the  effect 
of  the  superimposed  noise  component,  and  the  other  is  the 
different  selection  of  N and  t for  dimension  estimation.  In 
the  following  section,  we  will  investigate  the  noise  effect 
(filtering  effect)  on  dimension  estimation. 

Effect  on  Dimension  Estimation 
Noise  contamination  has  been  an  inevitable  problem  in 
collecting  EEG  data.  It  is  important  to  understand  the  role 
of  added  noise  in  the  quantified  dynamical  properties  of  the 
deterministic  signal  (80) . With  regards  to  the  experimental 


correlation  exponent 
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Figure  2-11  CEM  curves  for  two  15-second  sleep  two  EEG 
segments  of  the  same  patient. 
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data  like  EEG,  what  will  be  the  effects  of  the  external 
noise  fluctuations?  The  most  straightforward  approach  is  to 
investigate  the  effect  of  linear  filtering  on  the  dimension 
estimation  of  the  EEG  signal.  Filtering  actually  enhances 
coherence,  hence,  is  expected  to  decrease  the  correlation 
dimension.  To  get  a more  realistic  picture  of  how  filtering 
creates  spatial  correlation,  we  will  also  evaluate  the 
filtering  effect  on  the  correlation  dimension  of  white  noise 
sequences.  The  effects  of  filtering  bandwidth  and  filter 
types  are  investigated  as  well.  Hence,  this  study  involves 
two  steps:  (1)  evaluating  the  correlation  dimension  of  an 
original  and  filtered  pseudorandom  sequence,  and  (2) 
quantifying  the  effects  of  linear  filtering,  types  of 
filter,  and  filter  bandwidth  on  the  correlation  dimension  of 
the  EEG  segment. 

Three  linear  lowpass  digital  filters  (linear  phase  FIR, 
Chebyshev,  and  Butterworth)  are  designed  and  applied  to  the 
pseudorandom  sequence  and  EEG  signal.  Figure  2-12  to  2-14 
display  the  frequency  response  and  impulse  response  of  FIR 
(order  = 11,  first  zero  at  0.567T),  Chebyshev  (order  = 9, 
cutoff  frequency  = 0.4rr,  passband  ripple  = 0.04),  and 
Butterworth  (order  = 9,  cutoff  frequency  = 0.47T)  lowpass 
filters,  respectively.  A pseudorandom  number  sequence  with 
period  of  8192  samples  is  plotted  in  Figure  2-15  (the  top 
trace)  for  2500  samples,  which  are  generated  as  follows 
(81)  : 


magnitude  of  h(n) 


66 


0.25 
0.20 
0.15 
0.10 
0.06 
0.00 

Figure  2-12  FIR  lowpass  filter  (N=ll,  1st  zero  at  70  Hz) . 

(a)  Frequency  response;  (b)  Impulse  response. 
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Figure  2-13  Chebyshev  lowpass  filter  (N=9,  fc=50  Hz). 

(a)  Frequency  response;  (b)  Impulse  response. 
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Figure  2-14  Butterworth  lowpass  filter  (N=9,  fc=50  Hz).  • 
(a)  Frequency  response;  (b)  Impulse  response. 
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In+i  = (J  x In  + 1)  mod  M,  n = 0,  1,  M-l, 

(2.2) 

where  J = 4K  + 1,  M = 2L  (K  and  L are  integers) , and  initial 
condition  I0  = 1235.  The  sequence  is  then  filtered  by  a 
Butterworth  lowpass  filter  with  the  cutoff  frequency  of 
0.6477,  0.47 t,  and  0.1677  (the  other  three  traces  in  Figure  2- 
15) . Obviously,  the  random  sequence  is  much  smoother  (the 
bottom  trace  in  Figure  2-15)  after  filtered  by  a narrow-band 
(0.1677)  lowpass  filter,  which  shows  the  enhancement  of 
coherence.  In  the  following,  the  correlation  dimension  will 
be  estimated  for  the  four  sequences  to  verify  if  the 
correlation  estimation  detects  deterministic  components  in 
the  filtered  pseudorandom  noise. 

Filtering  Effect  on  the  Dimension  of  White  Noise 

The  MI  vs.  r curves  in  Figure  2-16  are  approximately 
flat  for  r > 2.  Theoretically,  an  infinite  white  random 
noise  sequence  is  completely  uncorrelated  and  has  infinite 
dimension.  Since  mutual  information  measures  the  general 
dependence  of  two  variables  (59) , the  MI  vs.  r curve  for  a 
random  sequence  should  be  zero  for  r f 0 and  1 for  r = 0. 

The  exceptionally  large  value  at  r = 1 is  due  to  the  finite 
period  (8192)  of  the  generated  pseudorandom  sequence.  The 
uncorrelated  property  of  a random  sequence  enables  us  to 
arbitrarily  choose  the  value  of  time  delay.  Here,  the  same 
parameters  (N  = 3750  and  r = 13)  will  be  used  for  both  the 
generated  random  sequence  and  EEG  segment. 


mutual  information 
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Figure  2-16  MI  curves  for  generated  white  noise. 
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The  CIM  curves  for  the  generated  pseudorandom  sequence 
are  plotted  as  solid  curves  in  Figure  2-17 (a)  to  2-17 (c), 
while  the  dashed  curves  represent  the  CIM  curves  for  three 
filtered  pseudorandom  sequences  (from  top  to  bottom,  fc  = 

80,  50,  and  20  Hz).  Figure  2-18  displays  the  CEM  curves  for 
a theoretical  white  noise  (straight  line)  and  the  four 
sequences.  Although  in  principle  no  correlation  dimension 
can  be  computed  due  to  the  nonsaturating  nature  of  the  CIM, 
Table  2-2  lists  the  correlation  exponents  estimated  for  n = 
16.  The  CIM  curves  for  the  original  pseudorandom  sequence 
manifest  no  parallelism  as  the  embedding  dimension 
increases,  which  is  reflected  in  the  nonsaturating 
characteristic  of  the  CEM  curve  (the  stars  and  solid  curve 
in  Figure  2-18) . Although  the  dashed  curves  in  Figure  2- 
17(a)  (fc  = 0.647T)  are  displaced  towards  smaller  ln(r)  , the 
slopes  are  basically  unchanged.  Accordingly,  the  CEM  curve 
(curve  3 in  Figure  2-18)  of  the  filtered  pseudorandom 
sequence  approximately  coincides  with  the  originally 
generated  pseudorandom  sequence  (curve  2 in  Figure  2-18) . 
Note  that,  for  each  specified  n,  the  solid  curve  of  Figure 
2-17  always  starts  off  with  a much  larger  ln(r)  than  the 
dashed  curve.  This  phenomenon  corroborates  the  dispersive 
nature  of  the  phase  trajectory  (Figure  2-19 (a)) 
reconstructed  from  random  noise,  which  does  not  converge 
towards  a subset  of  lower  dimensional  space.  Apparently, 
the  points  in  the  phase  trajectory  occupy  each  available 
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Figure  2-17  CIM  curves  for  original  and  filtered  white 
noise  (solid:  original,  dashed:  filtered), 
(a)  f c=80  Hz;  (b)  f c=50  Hz;  (c)  fc=20  Hz; 


correlation  exponent 
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Figure  2 18  CEM  curves  for  original  and  filtered  white 
noise. 
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Table  2-2 

Correlation  Exponent  for  Pseudorandom  Sequence 
(embedding  dimension  = 16) 


original  pseudorandom  seq. 

12.9007 

Butterworth 

(fc  = 80  Hz) 

11.6904 

Butterworth 

(fc  = 50  Hz) 

9.4885 

Butterworth 

(fc  = 20  Hz) 

7.7527 

Figure  2-19  Phase  portraits  constructed  from  a 3000-sample 
white  noise  with  delay  = 13. 

(a)  Original  white  noise; 

(b)  Filtered  white  noise  (fc=20  Hz) . 
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phase  space  dimension  and  resemble  a "random  walk"  in  phase 
space.  Therefore,  the  correlation  integral  ln[C(r)] 
vanishes  for  small  ln(r)  since  no  pair  of  data  points  can  be 
found  within  this  distance.  Besides  the  missing 
parallelism,  "knees"  no  more  occur  even  for  an  embedding 
dimension  as  high  as  16!  As  the  pseudorandom  sequence  is 
filtered  with  bandwidth  = 0.647T  and  0.47T  (the  2nd  and  3rd 
trace  in  Figure  2-15) , still  no  knee  is  explicitly  visible 
in  the  CIM  curves  (dashed  curves  in  Figure  2-17 (a)  and  2- 
17(b)).  However,  correlations  show  up  in  the  CEM  since  the 
correlation  exponent  does  not  increase  proportionally  to  the 
embedding  dimension  (Table  2-2) . When  more  severe  filtering 
is  applied  (bandwidth  = 0.16tt),  knees  show  up  in  the  CIM 
curves  (dashed  curves  in  Figure  2-17 (c) ) , which  indicates 
that  coherent  characteristics  due  to  the  filtering  begin  to 
appear  and  take  over  the  dynamical  activities.  Moreover, 
the  phase  trajectory  (Figure  2-19 (b) ) of  the  bottom  trace  in 
Figure  2-15  forms  an  attractor  with  finite  dimension,  which 
is  verified  by  the  tendency  of  saturation  of  the  CEM  curve 
for  n > io  (curve  5 in  Figure  2-18).  The  effect  of 
filtering  in  the  correlation  dimension  of  a white  noise  is 
evident  since,  with  embedding  dimension  16  (Table  2-2),  the 
correlation  exponent  decreases  from  11.69  (fc  = 80  Hz)  to 
7.75  (fc  = 20  Hz) . 

From  the  above  study  we  may  conclude  that  filtering,  a 
smoothing  operation,  creates  some  coherence  in  the  time 
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series  and  accordingly  the  correlation  dimension  is  reduced. 
Filtering  enhances  waveform  patterns  which  are  not  apparent 
in  the  original  signal.  This  side  effect  is  well  depicted 
in  Figure  2-15 (d)  and  2-17 (c)  (the  dashed  CIM  curves),  in 
which  case  a severe  filtering  induces  the  emergence  of  a 
subsystem  stated  before.  The  appearance  of  the  intermediate 
knee  in  CIM  curves  (dashed  curves  in  Figure  2-17 (c) ) is 
probably  caused  by  these  newly  observed  waves  and  patterns. 
Another  explanation  of  the  knee  phenomenon  in  the  severely 
filtered  random  sequence  might  be  due  to  the  clustering  of 
the  originally  evenly-distributed  amplitudes  into  small- 
amplitude  jittering  (average  = r^  and  large-amplitude 
oscillation  (average  = r2)  with  r2  <<  r2. 

Filtering  Effect  on  the  Dimension  of  EEG 

In  principle  the  correlation  dimension  estimates  the 
global  dynamical  property  of  a system,  therefore  it  is 
unlikely  to  be  sensitive  to  the  filter  types  which  may  cause 
a subtle  difference  in  the  waveform  patterns  of  a signal. 
While  some  approaches  of  signal  analysis,  such  as  waveform 
detection  and  recognition,  are  sensitive  to  the  local 
waveform  structure  produced  by  filters,  the  dimension 
estimate  is  expected  to  be  well  preserved  and  give  a 
consistent  result.  In  this  study,  we  are  investigating  the 
effect  of  filter  types  and  bandwidth  on  correlation 
dimension  of  the  EEG  signal.  The  sleep  stage  two  EEG  and 
its  filtered  versions  are  plotted  in  Figure  2-20  to  2-23. 
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Figure  2-20  30-second  sleep  stage  two  EEG. 
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Figure  2-21  30-second  sleep  stage  two  EEG  filtered  by 

FIR  lowpass  filter  (N=ll,  1st  zero  at  70  Hz).  1 sec 
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Figure  2-22  30-second  sleep  stage  two  EEG  filtered  by 
Chebyshev  lowpass  filter  (N=9,  fc=50  Hz). 
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Figure  2-23  30-second  sleep  stage  two  EEG  filtered  by 
Butterworth  lowpass  filter  (N=9,  f =50  Hz 
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High  frequency  components  with  smaller  amplitude  jittering 
are  significantly  filtered  out,  and  some  phase  shift  is 
introduced  into  the  filtered  signals.  The  segment  length  is 
set  at  15  seconds  (3750  samples) . To  show  how  filtering 
affects  the  structure  of  phase  trajectories,  Figure  2-24 
displays  the  2-D  projection  of  3-D  phase  trajectory 
constructed  from  the  3000-sample  segment  of  the  original  EEG 
(Figure  2-24 (a))  and  EEG  filtered  by  Butterworth  LPF  with  fc 
= 80  Hz  (Figure  2-24 (b) ) and  20  Hz  (Figure  2-24 (c) ) . 
Obviously,  the  zigzag  lines  of  Figure  2-24 (a)  are  smoothed 
out  so  that  a much  well-defined  structure  appears. 
Furthermore,  the  trajectory  curls  up  towards  the  center, 
which  results  in  the  contraction  of  the  high-density  region. 
In  other  words,  the  radius  of  the  basin  of  EEG  attractor 
seems  to  decrease. 

From  the  above  discussion,  it  appears  that  a reduction 
in  the  correlation  dimension  is  expected  for  a filtered  EEG 
signal.  However,  on  what  scale  does  the  correlation 
dimension  decrease? 

While  evaluating  the  filtering  effect  on  correlation 
dimension,  the  first  step  is  to  investigate  the  effect  on 
the  MI  curves  and  the  selection  of  time  delay.  Figure  2-25 
displays  the  MI  curves  for  the  original  EEG  segment  (stars  & 
solid),  the  segment  filtered  by  FIR  (stars  & dashed), 
Chebyshev  (circles  & dashed) , and  Butterworth  (diamonds  & 
dashed)  lowpass  filter  (the  first  15-second  segment  in 


84 


Figure  2-24 


v 

Phase  portraits  for  sleep  stage  two  EEG. 

(a)  Original  data;  (b)  fc=80  Hz;  (c)  fc=20  Hz. 
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Figure  2-25  MI  curves  for  filtered  EEG  segments. 
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Figure  2-20  to  2-23) . Obviously,  filtering  smooths  out  the 
MI  curve  and  removes  the  first  few  local  minima  (such  as  r = 
2,  6,  and  9 samples).  No  significant  local  minimum  is 
evidently  shown  in  the  MI  curves.  The  time  delay  is 
selected  according  to  the  rule  mentioned  before,  i.e.  the 
first  few  local  minima  having  the  magnitude  of  MI  on  the 
order  of  the  flat  value  of  MI  are  selected.  In  this  case 
the  magnitudes  of  MI  are  nearly  on  the  same  order  when  10  < 
r < 30.  Hence  time  delay  of  13  samples  is  chosen. 

In  Figure  2-26,  the  solid  curves  represent  the  CIM  for 
the  original  EEG  segment,  and  the  dashed  curves  for  the 
segment  filtered  by  FIR  (Figure  2-12),  Chebyshev  (Figure  2- 
13),  and  Butterworth  (Figure  2-14)  lowpass  filter.  The 
dashed  curves  are  distributed  over  the  range  of  smaller 
ln(r) , which  indicates  the  reduction  of  the  dynamical  extent 
as  observed  from  the  shrunk  attractors  (Figure  2-24).  Knees 
start  to  appear  at  a smaller  n,  which  is  n = 6 for  the  three 
filtered  EEG  segments  and  n = 14  for  the  unfiltered  one.  A 
decrease  of  correlation  dimension  by  2 is  observed  (Table  2- 
3).  As  shown  in  Figure  2-27,  quite  similar  estimations  are 
reached  for  the  three  filter  types.  Consequently,  the 
filter  type  does  not  seem  to  be  important  for  the  estimation 
of  correlation  dimension. 

To  investigate  the  effect  of  filter  bandwidth  on 
dimension,  Butterworth  lowpass  filter  with  9th  order  and  fc 
=80,  50,  and  20  Hz  is  applied  to  the  EEG  signal.  The 
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Table  2-3 

Correlation  Dimension  of  Original  and  Filtered  EEG 


signal 

correlation 

dimension 

original 

EEG 

6.559 

± 0.091 

filtered 

by  BW80a 

4.919 

± 0.066 

filtered 

by  BW50a 

4.706 

± 0.080 

filtered 

by  BW20a 

4.454 

± 0.066 

filtered 

by  CBS50b 

4.737 

± 0.075 

filtered 

by  FIR7  0C 

4.545 

± 0.065 

a)  BWn:  Butterworth  lowpass  filter  with  fc  = n Hz. 

b)  CBS50 : Chebyshev  lowpass  filter  (fc  = 50  Hz) . 

c)  FIR70:  FIR  lowpass  filter  (first  zero  at  70  Hz) . 
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Figure  2 


27  CEM  curves  for  original  and  filtered 
EEG  segments. 
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results  are  shown  in  Figure  2-28,  the  CEM  curves  for  all 
four  cases.  A reduction  in  bandwidth  does  cause  a slight 
decrease  of  the  correlation  dimension  (e.g.  a decrease  of 
0.21  when  bandwidth  is  narrowed  down  from  80  to  50  Hz). 
However,  the  difference  in  estimation  is  negligible  if  the 
filter  bandwidth  is  cautiously  selected,  i.e.  to  eliminate  a 
considerable  amount  of  noise  without  distorting  the 
deterministic  components. 

Discussion  and  Conclusion 

The  above  results  show  that  a reduction  in  bandwidth 
(from  80  to  20  Hz)  has  much  less  effect  on  correlation 
dimension  of  the  EEG  segment  than  the  one  of  random  noise 
seguence.  It  is  misleading  to  display  the  CEM  curve  of  the 
filtered  EEG  signal  in  comparison  with  the  theoretical 
straight  line  for  white  noise.  We  propose  that  the  CEM 
curve  of  the  EEG  signal  should  be  plotted  with  the  CEM  curve 
of  white  noise  for  the  same  measurement  bandwidth  (Figure  2- 
29)  . 

We  have  shown  that  the  severe  filtering  causes  the 
emergence  of  extraneous  patterns  in  the  random  noise 
sequence.  This  is  also  likely  to  happen  in  the  EEG  signal 
though  the  decrease  in  correlation  dimension  is  little 
(9.45%)  as  bandwidth  is  reduced  from  80  Hz  to  20  Hz.  As 
displayed  in  Figure  2-17 (c) , the  CIM  curves  of  a severely 
filtered  random  noise  seguence  exhibit  coherent 
characteristics.  This  phenomenon  is  also  observed  in  the 
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Figure  2-28  Effect  of  cutoff  frequency  on  CEM  curves 
for  sleep  stage  two  EEG. 


correlation  exponent 
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Figure  2-29 


CEM  curves  for  EEG  and  white  noise 

both  filtered  by  Butterworth  LPF  (fc=50  Hz) . 
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EEG  case — the  nonparallel  curves  in  the  lower  region  of  CIM 
curves  of  the  unfiltered  EEG  exhibit  parallelism  after 
filtered  with  a narrow  bandwidth  (20  Hz) . That  is  to  say, 
filtering  may  introduce  a subsystem  in  the  dynamics  since 
the  lower  region  of  CIM  characterizing  the  uncorrelated 
random  noise  has  been  converted  into  a region  showing 
coherence.  Conseguently , we  might  interpret  the  "knee" 
phenomenon  as  the  effect  of  two  subsystems  involved  in  the 
system  dynamics! 

From  the  above  experiment,  we  conclude  that  in  our 
study  of  the  EEG  nonlinear  dynamics  filtering  is  needed  to 
eliminate  the  effect  of  noise.  However,  further  research  is 
needed  to  investigate  the  effect  of  filtering  on  the 
dynamical  invariants.  From  the  CIM  curves  it  seems  that 
filtering  emphasizes  the  knee  phenomena,  displacing  it  not 
only  towards  smaller  r's  (as  expected),  but  also  towards 
lower  embedding  dimension.  This  last  point  may  indicate  the 
introduction  of  spurious  dynamical  features.  At  this  point 
of  the  research,  our  recommendation  is  to  cautiously  select 
the  filter  bandwidth  according  to  both  the  information  of 
the  signal  and  the  characteristics  of  CIM  curves,  so  that 
filtering  does  not  distort  the  system  dynamics. 

Dimensionality  Analysis  of  Lorenz  Model 
Why  do  we  want  to  analyze  the  dimensionality  of  the 
Lorenz  model  after  it  has  been  extensively  studied  for 
decades?  We  are  told  that  the  dimension  of  a Lorenz 
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attractor  with  system’s  parameters  of  a = 10.0,  r = 28.0, 
and  b = 8/3  is  2.05  ± 0.01  (17,64).  The  consistency  and 
accuracy  of  this  number  have  been  clearly  proven,  and  this 
is  the  importance  of  the  Lorenz  attractor:  it  can  be  used  as 
a benchmark,  to  test  new  experimental  procedures,  or  to 
answer  research  questions.  From  the  study  of  real 
biological  systems  we  realize  that  there  are  pre-requisites 
for  the  analysis  that  are  hardly  ever  met.  One  of  the 
lingering  questions  is  the  effect  on  the  quantified 
dynamical  properties  of  time  varying  system's  parameters. 

Are  we  capable  of  predicting  the  effect  of  time  varying 
system's  parameters  on  the  quantification  of  dynamical 
characteristics?  The  above  problems  cannot  be  neglected 
since,  in  our  study  of  brain  dynamics,  we  have  to 
realistically  assume  that  the  brain  is  a time-varying 
system. 

In  this  study,  the  dimension  estimate  of  the  Lorenz 
model  (41,52,54)  is  designed  to  illustrate  the  effect  of 
time  varying  system  parameters.  The  Lorenz  model  described 
by  the  following  equations 


will  be  numerically  simulated  by  Butcher's  fifth-order 
Runge-Kutta  method  (51) . Two  parameters  will  be  fixed  at  r 
= 45.92,  b = 4.0,  while  a varies.  With  the  initial 


dx/dt  = a (y  - x) , 
dy/dt  = x (r  - z)  - y, 
dz/dt  = xy  - bz, 


(2.3) 

(2.4) 

(2.5) 
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conditions  of  x(0)  = y(0)  = z(0)  = 1.0,  a transition  from 

stability  to  chaos  is  observed  when  a increases  above  6.3 

(one  decimal  digit  accuracy).  For  o between  6.3  and  39.8, 

the  Lorenz  model  keeps  on  exhibiting  chaotic  behavior. 

Above  this  value  it  becomes  stable,  i.e.  a fixed-point 

attractor,  and  so  forth  (Figure  2-30) . Accordingly,  the 

Lorenz  attractor  shows  intermittent  transient  chaos  (12) . 

The  study  involves  four  different  cases  of  parameter  (a) 

changing  during  the  evolution:  (1)  from  a stable  to  a stable 

system,  (2)  from  a stable  to  a chaotic  system,  (3)  from  a 

chaotic  to  a stable  system,  and  (4)  from  a chaotic  to  a 

chaotic  system.  One  factor,  the  rate  of  change  of  a,  plays 

an  important  role  in  the  experiment  since  whether  the  system 

dynamics  is  completed  or  not  affects  the  evaluation  of 

correlation  dimension.  Both  conditions,  including  the 

change  of  a before  and  after  the  settling  of  system 

dynamics,  will  be  investigated. 

Correlation  Dimension  of  Lorenz  Attractor  with 
Constant  Parameters 

Lorenz  attractors  for  o = 6,  10,  17,  24,  31,  and  40  (r 
= 45.92,  b = 4.0)  are  plotted  in  Figure  2-30(a)  to  2-30(f). 
The  pictures  show  the  top  view  (i.e.,  the  projection  of  the 
image  on  the  XY-plane)  of  a 3D-to-2D  parallel  projection, 
with  90°  rotation  angle  about  X-axis  and  0°  rotation  angle 
about  Y and  Z-axis.  The  stable  attractors  for  a = 6 and  40 
approach  the  same  fixed  point  (-13.41,  -13.41,  44.95)  at 
different  speed;  while  the  chaotic  Lorenz  attractors  for  a = 
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Figure  2-30  Lorenz  attractors  (r=45.92,  b=4.0,  a varies), 
(a)  a = 6;  (b)  a = 10;  (c)  a = 17; 

(d)  a = 24;  (e)  a = 31;  (f)  a = 40. 
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Figure  2-30  (continued) . 
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Figure  2-30  (continued) . 
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10,  17,  24,  and  31  track  a kind  of  double  spiral  in  three 
dimensions,  like  two  wings  of  a butterfly.  In  spite  of  the 
seemingly  infinite  complexity,  the  trajectories  always  stay 
within  certain  bounds  which  display  some  kind  of  order 
corroborating  the  chaotic  behavior.  However,  the  attractors 
with  different  cr's  have  distinct  boundaries  and  spatial 
extents.  The  dimensionality  therefore  will  reflect  how  an 
attractor  with  a specified  a spans  the  3D  space.  While 
evaluating  the  correlation  dimension  of  the  Lorenz  model,  we 
found  out  that  the  total  evolved  length  (in  seconds) 
required  for  obtaining  a stable  estimate  differs  especially 
between  a stable  and  chaotic  regime.  Another  considerable 
factor  to  be  investigated  is  the  step  size  (h)  used  in 
numerical  simulation. 

Theoretically,  a fixed-point  attractor  has  a dimension 
of  zero,  yet  our  study  of  the  stable  Lorenz  attractor  shows 
a non-zero  correlation  dimension  (Figure  2-31 (a)  and  2- 
31(f)).  Since  it  takes  a longer  time  (>  120  seconds)  to 
reach  and  finally  rest  on  the  fixed  point,  the  estimated 
correlation  dimension  converges  very  slowly.  On  the  other 
hand,  a chaotic  Lorenz  attractor  swirls  a large  area  in  the 
early  stages.  Even  with  a short  length  of  25  seconds,  we 
can  already  obtain  a convergent  estimation  of  /i  (Figure  2- 
31(b)  to  2-31 (e) ) . 

It  is  the  duration  in  seconds  that  actually  affects  the 
estimation.  To  save  the  computation  time,  we  estimate  the 
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Figure  2-31  (continued) 
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correlation  dimension  for  Lorenz  attractor  simulated  with 
different  step  sizes  (h) . As  a consequence,  we  find  that 
the  step  sizes  ranging  from  0.01  to  0.02  for  chaotic  Lorenz 
attractors  and  from  0.01  to  0.04  for  fixed-point  Lorenz 
attractors  are  appropriate  for  the  estimation  of  correlation 
dimension. 

Back  to  the  problem  of  "convergent  estimation  of  ju" , 
the  correlation  dimensions  of  Lorenz  attractors  for  o = 6.0, 

10.0,  17.0,  24.0,  31.0,  and  40.0  are  estimated  for  evolved 
length  =25,  50,  75,  and  100  seconds  (for  a = 6.0  and  40.0, 
length  = 120  seconds  is  also  analyzed).  Step  size  of  0.02 
(0.04  for  o = 6.0  and  40.0  with  120-second  duration)  is  used 
in  simulating  the  Lorenz  attractor.  Selecting  the  initial 
conditions  is  another  important  task  in  that  some  tiny 
differences  in  the  initial  state  could  lead  to  an  enormously 
different  outcome.  For  instance,  the  stable  Lorenz 
attractor  in  Figure  2-30 (f)  is  transformed  into  a chaotic 
attractor  (Figure  2-32)  if  the  system  originates  at  (1.0, 

1.0,  0.0).  in  this  study  the  attractors  start  from  an 
initial  position  (1.0,  1.0,  1.0).  The  ln[C(r)]  will  be 
computed  for  ln(r)  = o.l  to  3.4  (3.2  for  a = 6.0  and  10.0) 
stepped  by  0.1.  The  computation  time  required  in  this 
simulation  is  much  longer  for  a stable  Lorenz  attractor. 

For  example,  with  an  evolved  length  of  70.0  seconds  (N  = 

3500,  h = 0.02)  it  takes  218  minutes  for  a = 6.0  and  only  73 
minutes  for  o = 17.0,  i.e.  only  one  third  the  time  is  needed 
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Figure  2 32  Stable  Lorenz  attractor  becomes  chaotic 
with  initial  conditions  (1.0,1. 0,0.0). 
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when  compared  with  the  computation  of  the  stable  dynamics 
case!  The  Cl  curve  (Figure  2-31)  has  an  evident  linear 
Portion  for  ln(r)  from  0.1  to  2.7,  which  will  be  used  to 
compute  the  slope  (/x)  . The  estimated  correlation  dimension 
for  each  a with  four  evolved  lengths  is  shown  in  Table  2-4. 
The  listed  averages  and  standard  deviations  (in  %)  are  used 
to  verify  the  accuracy  (the  error  is  less  than  1%)  of  /x 
estimated  for  an  evolved  length  t > 50  seconds.  Both  the 
numerical  data  (/x  = 0.0978  and  0.1004  for  a = 6.0  and  10.0, 
respectively)  and  the  graphical  demonstration  (Figure  2- 
30(a)  and  2-30 (f))  show  that  the  Lorenz  model  for  a = 6.0 
approaches  a fixed  point  faster  than  the  one  for  a = 40.0. 
Apparently,  the  slope  of  the  Cl  continues  to  decrease  with 
increasing  t for  stable  attractors  (Figure  2-31 (a)  and  2- 
31(f)).  As  for  the  chaotic  attractors  (Figure  2-31(b)  to  2- 
31(e)),  the  overlap  of  Cl  curves  begins  for  t = 50  seconds 
and  above.  The  Lorenz  attractor  for  a = 17.0  (Figure  2- 
30(c))  spatially  extends  farther  than  the  other  chaotic 
regimes  and  conseguently  has  the  largest  correlation 
dimension  (2.0709).  When  a increases  above  17.0,  the 
dynamics  are  confined  to  a smaller  region  (Figure  2-30 (d) 
and  2-30 (e) ) , instead  of  tracing  orbits  of  varied  sizes. 

This  behavior  leads  to  a smaller  n (1.9299  for  a = 31.0)  for 
regimes  with  a > 17.0.  The  silhouette  of  an  attractor 
seemingly  reflects  the  estimated  correlation  dimension. 
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Table  2-4 

Correlation  Dimension  of  Lorenz  Attractors 
with  Constant  Parameters 


6.0 

10.0 

a 

17.0 

24.0 

31.0 

40.0 

25a 

0.8704 

2.1292 

1.9460 

2.0224 

1.9481 

1.0200 

( 15)  b 

(8) 

(9) 

(9) 

(9) 

(15) 

50 

0.2809 

2.0493 

2.0302 

1.9702 

1.9503 

0.2926 

(92) 

(33) 

(35) 

(33) 

(36) 

(82) 

75 

0.1682 

2.0448 

2.0607 

1.9841 

1.9357 

0.1711 

(218) 

(69) 

(73) 

(75) 

(76) 

(211) 

12  0cd 

0.0978 

2.0267 

2.0709 

1.9865 

1.9299 

0.1004 

100 

(271) 

(123) 

(129) 

(133) 

(136) 

(260) 

averagee 

2.0403 

2.0539 

1.9803 

1.9386 

standard 

deviation 

t 

0.479% 

0.842% 

0.363% 

0.443% 

a)  evolved  length  in  second. 

b)  The  number  in  the  parentheses  shows  the  computation  time 
in  minutes  for  estimating  Cl. 

c)  h = 0.03  and  N = 4000  for  a = 6.0  and  40.0. 
h = 0.02  and  N = t/h  for  other  a's. 

d)  evolved  length  t = 120  seconds  for  a=6.0  and  40.0, 
and  t = 100  seconds  for  other  a's. 

e)  The  correlation  dimension  for  t = 25  is  excluded  in  the 
calculation  of  average. 

f)  This  item  represents  the  standard  deviation  as  the 
percentage  of  the  average. 
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Correlation  Dimension  of  Lorenz  Attractor  with 
Changing  Parameters 

As  described  in  the  beginning,  four  cases  of  parameter 
(ct)  variation  will  be  studied.  The  initial  position  is 
(1.0,  1.0,  1.0).  We  investigate  the  simplest  case,  i.e.  a 
is  assumed  to  change  once  during  the  evolution.  The  a of  a 
stable  (chaotic)  Lorenz  attractor  will  be  called  a stable  a 
(a  chaotic  a) . The  attractor  of  a Lorenz  model  with  a 
varying  from  m (the  first  t1  seconds)  to  n (the  final  t2 
seconds)  will  be  called  an  m-*n  (t1,t2)  attractor.  Figure  2- 
33  shows  examples  of  Lorenz  attractors  with  time  varying  a, 
the  31^40  (35,35)  and  31^17  (35,35)  attractor. 

A step  size  of  0.02  seconds  is  used  to  simulate  the 
ct1-*ct2  Lorenz  model  under  all  regimes.  Table  2-5  lists  the 
estimated  correlation  dimension  for  o1-*a2  (35,35)  attractors 
(excluding  those  with  superscript) , and  the  Cl  curves  are 
plotted  in  Figure  2-34.  The  simplest  case  in  this  study  is 
the  system  with  a parameter  changing  between  stable  values. 
The  6->40  and  40->6  (35,35)  attractors  look  like  those  with 
constant-a  systems  (Figure  2-30(a)  and  2-30(f)).  Hardly  any 
difference  in  correlation  dimension  (Figure  2-34 (a))  is 
observed. 

It  seems  that  the  Lorenz  model  with  a changing  from  a 
chaotic  to  a stable  (chaotic)  number  is  stable  (chaotic) . 
However,  it  is  quite  another  story  after  further 
investigation.  How  an  m-*n  (t^,t2)  attractor  behaves 
actually  depends  on  t-^  and  t2.  For  instance,  a 40-+31 
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Figure  2-33  Lorenz  attractors  with  a varying  (r=45.92 

(a)  The  31-*40  (35,3  5)  attractor; 

(b)  The  31-*T7  (35,35)  attractor. 


109 


Table  2-5 

Correlation  Dimension  of  Lorenz  Attractor 
with  a Varying 


to 

6.0 

a varies  from 
10.0  17.0 

31.0 

40.0 

6.0 

1.6208 
( c-*-s ) a 

0.1866 

(S^S) 

10.0 

1 . 3883b 
0.4692 

(S^C) 

2 . 0556c 
2.2371 
(C-+C) 

17.0 

2.2919c 

2.2450 

(C^C) 

2.1622° 

2.1275 

(C^C) 

1 . 8016d 
0.5003 
(s-C) 

31.0 

2.0405° 

2.0405 

2.1182e 

(C^C) 

1.1513° 

0.1865 

(S^C) 

40.0 

0.1829 

(S^S) 

1 . 0436f 
0.6829 
(C-S) 

a) 

b) 

c) 

d) 

e) 

f) 


c:  chaotic,  s:  stable. 

6->l0  (20.6,49.4)  attractor. 
al*°2  (17.55,17.55)  attractor. 
40-*17  (20,50)  attractor. 
i7->31  (50,50)  attractor. 

31-*40  (20,50)  attractor. 


ln(C(r) ) 
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Figure  2-34  Cl  curves  for  Lorenz  attractors  with  a varying 
(r=45 . 92 , b=4 . 0 , h=0. 02) . 

(a)  Stable  a -*  stable  a attractors; 

(b)  Chaotic  a -*  chaotic  a attractors; 

(c)  Chaotic  a -*■  chaotic  o attractors; 

(d)  Stable  a * — ♦ chaotic  a attractors; 

(e)  Stable  a * — * chaotic  a attractors. 


ln(C(r))  s ln(C(r)) 
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Figure  2-34 


(continued) 


ln(C(r))  s ln(C(r) ) 
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Figure  2-34  (continued) 
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attractor  is  a fixed  point  like  the  one  shown  in  Figure  2- 
30(f)  if  the  variation  of  a occurs  after  an  evolved  length 
tx  = 35  seconds;  however,  the  system  becomes  chaotic  if 
variation  occurs  earlier,  for  instance  at  tj_  = 16  seconds. 
This  phenomenon  is  due  to  the  intrinsic  dependency  on  the 
initial  conditions  of  a chaotic  system.  Since  the  final 
state  of  the  c^-system  (a  Lorenz  model  with  a = a2)  is  the 
initial  position  of  the  02-system,  t1  biases  the  ultimate 
behavior  of  a cr1-+a2  attractor.  Hence  for  a Lorenz  model 
with  o varying  from  a chaotic  to  a stable  a,  the  attractor 
does  not  necessarily  converge  toward  a fixed  point,  and  vice 
versa. 

The  third  case  examined  is  the  dimension  of  an 
attractor  with  chaotic  o1  and  a2  (Figure  2-34 (b)  and  2- 
34(c)).  For  chaotic  Lorenz  attractors,  an  evolved  length  of 
35  seconds  is  enough  for  a convergent  estimation  of  /x. 
Therefore,  a 0^02  (35,35)  attractor  with  chaotic  ctx  and  ct2 
displays  an  interweaved  image  of  both  the  and  02- 
attractor  (if  the  final  state  of  the  a^-attractor , after 
evolving  for  35  seconds,  will  not  drive  the  02-attractor  to 
a fixed  point) . The  correlation  dimension  increases  in 
comparison  with  those  for  constant-a  attractors  due  to  the 
increase  of  randomness  of  the  model  with  parameter 
variation.  The  numbers  in  parentheses  listed  below  indicate 
the  percentages  of  increase  (or  decrease,  if  a negative 
number)  of  the  dimension  of  a Q\*g2  system  with  respect  to  a 
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ax  and  CT2-system. 

10-17 17-10 17-31 31-17 

% (10.77,8.41)  (8.03,10.38)  (-1.47,5.73)  (10.24,2.73) 

The  negative  number,  which  manifests  the  decrease  of  /z 
(2.0405)  of  a 17—31  (35,35)  attractor  when  compared  to  the  /z 
of  the  attractor  with  a = 17  (2.0709),  calls  our  attention 
to  further  investigation.  With  100  seconds  of  evolved 
duration,  however  a larger  /z  (2.1182)  is  obtained  for  a 
17—31  (50,50)  attractor.  While  considering  a system  with 
the  parameter  varying  from  a chaotic  to  a stable  (or  a 
stable  to  a chaotic)  a,  diverse  values  of  n (Figure  2-34 (d) 
and  2-34 (e) ) emerge  due  to  the  complicated  combinations. 

The  phase  trajectory  could  be  a fixed-point  attractor  having 
an  almost  zero  dimension,  such  as  the  40-31  (35,35) 
attractor  (/z  = 0.1865),  or  a butterfly-like  chaotic 
attractor  (the  40-17  (20,50)  attractor,  n = 1.8016). 
Discussion  and  Conclusion 

After  analyzing  the  effects  of  a time-varying  parameter 
on  the  correlation  dimension,  the  problem  of  whether  the 
parameter  variation  produces  a large  increase  in  the 
dimension  (if  /z  does  increase)  is  unsolved.  Even  for  a 
17"+31->17  (24,24,22)  system  with  a varying  twice,  the 
correlation  dimension  (2.0971)  shows  only  a 1.27%  increase. 

A higher  rate  of  increase  in  n may  occur  in  such  case  as  a 
10-17-31  (24,24,22)  system  (/z  = 2.2688),  although  we  are  not 
aware  of  how  or  in  what  manner  the  system's  parameter 
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changes  by  studying  the  dimension.  Yet,  from  this  study  we 
notice  that  a complicated  system  with  a parameter  changing 
indefinitely  . . .-*an-*  •••)  is  expected  to  rise  its 

dimension  close  to  the  system's  degrees  of  freedom  (i.e., 
the  embedding  dimension).  For  instance,  a 17-*-31-*17-*-31-»  ... 
(0.4, 0.4, 0.4, 0.4,  ...)  system,  i.e.  periodically  (period  = 
0.4  second)  varying  the  a between  two  values  17  and  31  for 
70  seconds,  has  an  estimated  correlation  dimension  of  2.1602 
(4.32%  increase).  And  the  model  with  a increasing  by  a step 
of  0.2  every  0.525  second  shows  a 12.04%  increase  in  n 
(2.2961).  More  complicated  pattern  of  parameter  variation 
may  result  in  an  attractor  having  the  correlation  dimension 
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CHAPTER  3 

ESTIMATION  OF  THE  LARGEST  LYAPUNOV  EXPONENT 
Experimental  Implementation  and  Considerations 
It  has  been  shown  (25,30,69)  that  the  spectrum  of 
Lyapunov  exponents  provides  the  most  useful  tool  for 
examining  the  dynamical  properties  of  chaotic  systems.  The 
fine  detail  of  the  attractor  can  be  revealed  through  the 
Lyapunov  exponents,  therefore  this  technigue  has  the 
potential  to  discriminate  between  deterministic  chaotic 
behavior  and  pure  random  activity.  Wolf  et  al.  (16) 
developed  a technigue  which,  for  the  first  time,  estimated 
the  two  largest  non-negative  Lyapunov  exponents  from  finite 
experimental  data  segments.  There  are  other  methods 
developed  to  calculate  Lyapunov  exponents,  however,  only  for 
systems  with  known  eguations  of  motion  (29,30,82).  The 
basic  concept  of  the  technique  developed  by  Wolf  et  al. 
(16,70)  originated  from  the  works  of  Shimada  and  Nagashima 
(29)  and  Bennetin  et  al.  (30)  . It  will  be  applied  here  to 
provide  more  knowledge  of  the  structure  of  an  EEG  chaotic 
attractor.  As  far  as  we  know  only  one  published  result  of 
EEG  Lyapunov  exponent  exists  (8) . However,  a reliable  and 
satisfactory  result  is  difficult  to  reach  due  to  some 
inevitable  problems  in  experimental  data  and  too  many 
undetermined  parameters  when  implementing  the  technique. 
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First  we  will  bring  forth  some  problems  encountered  in 
the  estimation  of  Lyapunov  exponents.  To  quantify  chaotic 
and  deterministic  behavior  in  experimental  data,  it  is 
necessary  to  have  a reasonable  quantity  of  "accurate"  data. 
Especially  for  the  estimation  of  Lyapunov  exponents  which  is 
highly  sensitive  to  even  a minute  variation  of  the 
attractor's  structure,  the  quality  of  experimental  data  is 
considerably  important.  Although  an  enormous  amount  of  the 
EEG  data  can  be  collected  from  the  continuous,  long-term 
recording,  the  accuracy  or  consistency  of  the  data  is  hardly 
validated  for  the  complex  brain  activities.  The  additive 
random  noise  and  the  time-varying  nature  of  the  EEG  dynamics 
are  the  most  evident  factors  which  affect  the  reliability  of 
exponent  estimation.  The  estimation  may  be  misleading 
without  recognizing  the  change  of  brain  state.  Moreover, 
some  uncertain  factors  in  the  technique  still  persist, 
further  complicating  the  estimation. 

The  implementation  parameters  required  include: 
duration  of  signal  (N  in  samples) , embedding  dimension  of 
the  reconstructed  trajectory  (n) , time  delay  for  space 
reconstruction  (t  in  samples) , sampling  rate  (fs  in  Hz) , 
maximum  length  scale  limitation  (SCALMX  or  pscalmx) , minimum 
length  scale  limitation  (SCALMN  or  pscalmn) , and  a constant 
evolved  step  between  replacement  attempts  (EVOLV  in 
samples) . The  major  task  of  the  algorithm  is  to  carry  out 
repeated  cycles  of  propagating  and  replacing  the  principal 
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axis  vector,  (P(t),Q(t'))  and  (P (t+EVOLV) ,Q (t") ) in  Figure 
3-1.  Here  we  introduce  the  term  "local  structure"  for 
convenience.  As  illustrated  in  Figure  3-1,  the  local 
structure  depicts  the  relative  geometric  connection  between 
a segmental  fiducial  trajectory  (i.e.  the  segment  from  P(t) 
to  P (t+EVOLV)  on  the  fiducial  trajectory)  and  another 
trajectory  with  an  initial  position  Q(t')  that  is  the 
spatial  nearest  neighbor  of  the  fiducial  point  P(t)  and  is 
also  evolving  for  EVOLV  samples  (Q (t ' +EVOLV) ) . The  largest 
Lyapunov  exponent  T!  is  estimated  (16,83)  by  averaging  the 
divergent  or  convergent  rates  of  all  the  local  structures 
traversing  an  entire  attractor.  The  divergent  or  convergent 
rate  of  one  local  structure  (will  be  called  "instantaneous 
running  exponent,"  irexp)  is  calculated  as  follows 
(Appendix) : 

irexp  = ln(df/di)  / [EV0LV*dt*ln(2) ] , (3.1) 

where  dt  is  the  sampling  interval  in  second  (dt  = l/fs) , and 
di  and  df  are  the  initial  and  final  separation, 

D[P(t) ,Q(t*) ] and  D[P (t+EVOLV) ,Q(t'+EVOLV) ] . Initial 
separation  di  of  any  local  structure  cannot  be  smaller  than 
SCALMN  that  defines  the  length  scale  on  which  noise  is 
expected  to  appear.  If  final  separation  df  is  larger  than 
SCALMX  (an  estimate  of  the  length  scale  on  which  the  local 
structure  of  the  attractor  is  no  longer  probed) , a 
replacement  (Q(t")  in  Figure  3-1)  for  Q(t'+EVOLV)  is  needed. 
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Figure  3-1  Local  structure  defined  in  exponent  estimates. 
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Otherwise,  (P (t+EVOLV) , Q (t ' +EVOLV) ) will  become  the 
principal  axis  vector  of  the  next  local  structure,  and  df  is 
its  initial  separation  di.  When  a replacement  is  needed, 
the  search  for  a best  point  (e.g.  Q(t"))  will  be  guided  by 
two  conditions:  (1)  SCALMN  < separation  D[P(t+EVOLV) ,Q(t") ] 

< SCALMX , and  (2)  very  importantly,  the  orientation  change 
must  be  as  small  as  possible.  Here  the  "orientation  change" 
indicates  the  angle  between  the  old  vector 
(P (t+EVOLV) ,Q(t'+EV0LV) ) and  the  replacement  vector 
(P (t+EVOLV) ,Q(t") ) . The  detailed  documentation  and  programs 
are  available  in  an  internal  report  of  the  Computational 
Neuroengineering  Laboratory,  University  of  Florida. 

In  the  following  we  discuss  the  data  conditions  and  the 
implementation  parameters  used  in  the  algorithm  (Appendix) 
evaluating  r^.  The  sleep  stage  two  EEG  (sampling  rate  = 250 
Hz)  is  still  our  object  for  the  estimation  of  the  largest 
Lyapunov  exponent  (r^)  . 

Embedding  Dimension  n and  Time  Delay  r 

These  two  parameters  can  be  determined  beforehand  by 
the  methodology  of  dimension  estimate  introduced  in  Chapter 
2 • From  the  MI  vs . t and  CIM  curves , an  adequate  t and 
number  of  embedding  variables  n can  be  determined. 

Duration  of  EEG  Signal 

The  correlation  dimension  of  sleep  stage  two  EEG  begins 
to  converge  for  N > 3750,  i.e.,  the  estimation  inaccuracy  is 
tolerable  if  a segment  length  of  no  less  than  3750  samples 
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is  used.  Generally,  the  correlation  dimension  is  less 
sensitive  to  the  little  variation  of  local  structures  of  the 
phase  trajectory  since  it  quantifies  the  global  property  of 
a system  (distance  between  pairs  of  points  irrespective  of 
ordering) . However,  the  Lyapunov  exponents  measure  the 
exponential  rates  of  divergence  or  convergence  of  nearby 
orbits  and  accordingly  is  much  affected  by  the  subtle  local 
orbit  structure  at  any  moment.  In  addition,  the  algorithm 
(Appendix)  operates  in  the  manner  that,  after  each  evolution 
step  (by  EVOLV  samples)  of  the  fiducial  trajectory,  a "local 
structure"  is  established  by  choosing  a nearest  point  to  the 
ongoing  segmental  fiducial  trajectory.  The  nearest  point  is 
chosen  from  all  the  available  data  points  used  to  construct 
the  attractor . Hence  the  amount  of  data  points  may  bias  the 
formation  of  local  structures  and,  consequently,  the 
estimation  of  in  this  study,  a segment  length  of  15 

seconds  is  first  used  as  a standard  for  sleep  stage  two  EEG, 
then  the  result  will  be  compared  with  those  evaluated  from 
different  durations  (4,  8,  30,  45,  and  60  seconds).  At  the 
same  brain  state,  different  segments  are  analyzed  to  test 
the  consistency  and  reliability  of  the  estimate. 

Evolution  Time  EVOLV  between  Replacement 

From  Figure  3-1,  a segmental  fiducial  trajectory 
represents  the  evolution  of  fiducial  point  from  P(t)  to 
P ( t+EVOLV) . Note  that  the  t and  EVOLV  are  measured  in  terms 
of  multiples  of  data  "sample"  since  we  are  dealing  with  a 
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discrete  signal.  Picking  a large  EVOLV  considerably  reduces 
the  computation  time  needed.  Nevertheless,  it  may  make  the 
volume  element  of  the  local  structure  (P,Q)  grow  excessively 
big.  In  addition,  too  large  an  EVOLV  will  cause  an 
erroneous  estimation  of  since  the  two  trajectories  might 
pass  through  a folding  region.  On  the  other  hand,  much 
orientation  error  will  be  introduced  into  the  estimation  if 
the  EVOLV  is  too  small  (16).  Here  the  "orientation  error" 
indicates  the  accumulated  angular  changes  resulted  from  the 
replacements  of  Q(t'+EVOLV).  Wolf  et  al.  have  proposed  that 
an  EVOLV  in  the  range  0.5  to  1.5  orbits  almost  always 
provides  stable  estimation  of  the  Lyapunov  exponents. 

However,  in  the  case  of  EEG's,  the  mechanism  of  system 
dynamics  is  unknown,  not  to  mention  the  orbital  size! 
Estimation  of  the  orbital  size  of  the  EEG  attractor 

To  picture  how  the  signal  itself  can  be  associated  with 
the  orbital  activities  of  a phase  trajectory,  we  propose  to 
the  information  from  the  power  spectrum  of  the 
signal  since  the  regularities  of  the  time  series,  to  our 
belief,  sway  the  revolving  motion  of  the  attractor.  We 
present  a short  experiment  below  to  explain  our  model. 

First  of  all,  the  power  spectrum  of  sleep  stage  two  EEG 
(4096  samples)  filtered  by  Butterworth  lowpass  filter  (fc  = 

30  Hz)  is  calculated  with  the  freguency  resolution  of 
250/1024  (i.e.,  0.24414)  Hz.  For  simplification,  we  will 
not  use  all  the  512  sinusoids  (512  frequencies)  to  simulate 
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the  EEG  segment;  instead,  only  those  frequency  components 
with  power  larger  than  a specified  threshold  (e.g.  50%,  25%, 
and  10%  of  the  maximum  power)  are  selected.  The  discrete 
time  series  is  then  represented  in  terms  of  sinusoidal 
sequences  involving  those  frequency  components.  The 
results,  three  1250-sample  segments,  are  plotted  in  Figure 
3-2  with  the  power (frequency)  information  listed  in  Table  3- 
1.  The  three  segments  in  Figure  3-2  will  be  called  segment 
A,  B,  and  C respectively.  Apparently,  the  original  signal's 
power  is  mostly  distributed  in  the  range  of  frequency  less 
than  4.0  Hz.  Although  segment  B involves  more  than  twice 
the  frequency  components  of  segment  A,  we  observe  the 
similarity  between  these  two  segments,  in  which  the 
waveforms  oscillate  less  than  4 cycles  each  second.  The 
mutual  information  of  the  three  segments  reveals  an 
interesting,  yet  predictable,  phenomenon.  The  MI  curve  of  a 
signal  composed  of  higher  frequencies  (segment  C)  descends 
more  steeply  to  the  first  local  minimum,  which  occurs  much 
ear-^er  (a*-  a smaller  delay)  than  the  one  composed  of  only 
lower  frequencies  (segment  A) . This  is  expected  since, 
between  the  original  and  delayed  time  series,  the  degree  of 
correlation  or  dependency  varies  rapidly  for  a high- 
frequency  signal,  and  consequently  the  delay  coordinates 
should  be  constructed  with  a small  r.  For  segment  A and  B, 
the  values  for  delay  are  observed  to  be  the  same  (18 
samples) , which  is  12  for  segment  C.  The  attractors  are 


125 


Table  3-1 

EEG  Segment  Simulation  by  Sinusoids 


Freauencies  (Hz) 

magnitude  of 
50%  25% 

power 

10% 

0.0000 

— 

— 

1.1842 

0.4883 

— 

— 

1.1469 

0.7324 

— 

— 

2.2924 

0.9766 

— 

3.3929 

3.3929 

1.2207 

— 

4.7764 

4.7764 

1.4648 

8.3697 

8.3697 

8.3697 

1.7090 

— 

3.4472 

3.4472 

1.9531 

7.6621 

7.6621 

7.6621 

2.1973 

7.2583 

7.2583 

7.2583 

2.4414 

10.0000 

10.0000 

10.0000 

2.6855 

— 

4.2155 

4.2155 

2.9297 

6.4895 

6.4895 

6.4895 

3 . 1738 

5.8738 

5.8738 

5.8738 

3.4180 

— 

4.5589 

4 . 5589 

3.6621 

— 

2.6343 

2 . 6343 

3.9063 

— 

3.9039 

3 . 9039 

4.1504 

— 

1.0828 

4.3945 

— 

_ 

2 .2100 

4.8828 

— 

— 

1.6699 

5.1270 

— 

_ 

2 . 3988 

5.3711 

— 

_ 

2 .4036 

5.8594 

— 

_ 

1.5207 

6.3477 

— 

_ 

1.3455 

6.8359 

— 

_ 

1.2275 

8.7891 

— 

_ 

1.0607 

10.9863 

— 

_ 

1.0653 

13.9160 

— 

_ 

1 . 0465 

14.1602 

— 

1 . 5644 

14.4043 

— 

1. 1429 

14.6484 

- 

- 

1.0686 
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reconstructed  in  a 3-dimensional  phase  space  and  projected 
in  the  plane.  Figure  3-3  only  shows  the  trajectories 
evolving  for  1 second  (the  range  from  0.2  to  1.2  second  in 
Figure  3-2) . The  attractor  A and  B basically  have  the  same 
structure  and  traverse  the  space  for  the  same  number  of 
orbits  of  different  sizes.  However,  for  only  a short 
duration  the  trajectory  of  segment  C already  gets  into  a 
mesh  of  unresolved  structure.  The  small-amplitude  up-and- 
down  ripple  in  the  time  series  contributes  to  the  small 
orbits  having  the  fast  and  frequent  change  of  direction 
along  the  corresponding  axis.  On  the  other  hand,  the  motion 
about  the  axis  constructed  from  the  slowly  varying  waveforms 
changes  smoothly  and  slowly. 

The  example  of  a 3-D  attractor  described  thus  far 
seemingly  oversimplifies  the  problem  encountered. 
Nevertheless,  the  concept  is  the  same  in  principle  for  our 
complicated  EEG  attractor  which  involves  more  than  five 
degrees  of  freedom.  It  is  obvious  that  the  occurrence  of  a 
peak  in  the  time  series  indicates  a change  in  direction 
along  one  of  the  axis  in  phase  space.  Thinking  of  an  n-D 
attractor  reconstructed  from  the  original  signal  and  (n-1) 's 
delayed  versions,  the  trajectory  traverses  a space  described 
by  the  coordinate  (alfa2, . . . ,an)  . At  time  tp,  a peak 
appears  on  the  ith  time  series,  i.e.,  ( i— 1 ) th  delayed 
version,  and  as  a consequence  the  trajectory  is  observed  to 
have  a directional  turn  along  a^axis  at  that  moment,  which 
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(c) 


Figure  3-3 


Trajectories  constructed  from  segments  in 
Figure  3-2  with  r = 18  (a) , 18  (b) , and  12 


(c)  . 
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may  cause  a folding  or  stretching. 

From  our  illustration  up  to  this  point,  we  still  have 
no  criterion  for  determining  the  orbital  size.  However  from 
the  'peaks'  information,  we  may  figure  out  how  frequently 
should  the  replacement  in  the  algorithm  be  executed  to  avoid 
the  folding  process.  In  addition,  it  is  worth  investigating 
how  the  "place"  of  replacement  affects  the  estimated 
exponent.  The  Lorenz  model  is  once  again  used  to  exemplify 
the  procedure.  The  1^  of  the  Lorenz  attractor  (o  = 16.0,  r 
= 45.92,  b = 4.0)  (54)  plotted  in  Figure  3-4 (a)  is  estimated 

under  three  variable-EVOLV  conditions  concerning  the 
replacements  which  (1)  always  take  place  at  the  peaks  of  the 
time  series  X(t),  Y(t),  and  Z(t)  (Figure  3-4 (b) , 3-4 (c),  and 
3—4 (d) ) , (2)  avoid  all  the  peaks,  and  (3)  occur  partially  at 
the  peaks.  Butcher's  fifth-order  Runge-Kutta  method  (51)  is 
applied  to  simulate  the  phase  trajectory  with  a step  size  of 
0.01  second,  N = 3750,  and  an  initial  condition  X(0)  = Y(0) 

= -15.0,  and  Z (0)  = 47.0.  The  largest  Lyapunov  exponent  is 
2.81,  1.37,  and  1.73  (in  (16),  r1  = 2.16  bits/s) 
respectively  for  the  three  cases  described  above.  Our  study 
shows  that  there  may  be  lots  of  orientation  change  of  the 
phase  trajectory  involved  in  case  (1),  which  leads  to  an 
overestimate  of  r^*  on  the  other  hand  (case  (2)),  the  rx  is 
underestimated  since  the  small  scale  local  structure  have 
involved  the  folding/stretching  region.  Therefore,  we 
decided  to  use  a constant  EVOLV  which  will  contemplate  both 
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Figure  3-4  Lorenz  attractor  (r=45.92,  b=4.0,  a=16.0). 

(a)  2-D  projection  of  3-D  trajectory; 

(b)  X-axis  time  series;  (c)  Y-axis  time  series; 
(d)  Z-axis  time  series. 
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Figure  3-4  (continued) . 
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replacement  strategies. 

Length  Scales  Limitation  (SCALMX  and  SCALMN) 

It  is  necessary  to  impose  a limit  on  the  scales  used  to 
define  each  evolved  local  structure  (P,Q) . SCALMN  and 
SCALMX  place  restrictions  on  the  distance  between  the 
fiducial  point  P(t)  and  the  selected  nearest  point  Q(t'). 
SCALMN  defines  the  length  scale  on  which  noise  is  expected 
to  appear;  and  SCALMX  confines  the  length  scale  of  the  local 
structure  being  investigated.  In  addition,  SCALMX  sets  the 
maximum  allowable  distance  at  the  end  of  each  EVOLV-sample 
propagation,  where  a replacement  (of  Q(t'+EVOLV))  is 
attempted  when  the  final  separation  of  the  evolved  pair  is 
greater  than  SCALMX.  Since  those  points  closer  than  SCALMN 
will  be  rejected  at  each  replacement  time,  the  effect  of 
noise  on  exponent  estimates  may  be  further  reduced.  SCALMN 
and  SCALMX  are  derived  from  the  following  equations 
(Appendix) : 

SCALMN  = pscalmn  x extent,  (3.2) 

SCALMX  = pscalmx  x extent,  (3.3) 

where  extent  is  the  largest  interpoint  distance  in  the  n- 
dimensional  attractor,  and  pscalmn  (pscalmx)  specifies  the 
percentage  of  extent  by  which  the  SCALMN  (SCALMX)  is 
defined.  The  pscalmn  (pscalmx)  from  1%  to  15%  (10%  to  40%) 
will  be  utilized  to  investigate  the  possible  range  of  stable 
exponent  estimates  and  the  effects  on  the  estimation  of  r^. 
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Experimental  Results 

The  dimension  and  time  delay  represent  the  essential 
knowledge  for  reconstructing  the  attractor  from  a time 
series.  From  the  MI  curves  of  six  15-second  segments,  the 
time  delay  of  9 samples  is  selected  as  before.  From  the 
dimensional  analysis,  we  are  able  to  determine  the  minimum 
number  of  variables  n necessary  for  describing  the  dynamics 
of  the  EEG  segments.  The  correlation  dimensions  of  six 
consecutive  15-second  segments  (7.237,  6.282,  7.051,  7.103, 
7.054,  and  7.317)  show  a slight  deviation  (4.8%)  which  may 
result  from  the  variation  of  system's  parameters  (discussed 
in  Chapter  2 and  Chapter  5) . To  our  belief,  a 90-second 
segment  categorized  as  a stage  two  sleep  EEG  is  very  likely 
to  involve  several  sets  of  system  parameters  that  may  not  be 
critical  for  the  EEG  visual  analysis,  but  that  can  cause 
different  estimations  in  highly  sensitive  methods  like  the 
Lyapunov  exponents.  Furthermore,  the  values  of  embedding 
dimension  at  which  parallelism  of  the  CIM  curves  develops 
are  not  all  the  same  (n  = 10,  9,  9,  10,  9,  and  9 for  these 
six  consecutive  segments) , and  for  each  case  the  appropriate 
n will  be  used.  Note  that  the  time  series  used  here  is 
noisy,  which  may  cause  undesired  effects  in  this  study  since 
the  estimation  of  Lyapunov  exponents  is  very  sensitive  to 
noise.  Moreover,  the  implementing  parameters  may  vary  with 
the  signal's  conditions  (i.e.  the  raw  signal  and  the 
filtered  signal) . 


135 


Effect  of  EEG  Segment  Selection  in  Estimation 

The  evolved  step  EVOLV  = 0.06  second  (15  samples)  is 
used  in  this  study  with  the  scale  limitations  pscalmx  =20% 
and  pscalmn  = 5%.  Figure  3-5  plots  the  running  exponents 
vs.  evolution  time  curves  (which  will  be  called  "Lexp-t 
curves”  hereinafter)  for  six  consecutive  15-second  segments 
(total  90-second  duration) . "Running  exponent"  (Lexp) 
represents  the  averaged  divergent/ convergent  rate  at  any 
moment  during  the  evolution  of  the  fiducial  trajectory 
(e.g.,  the  running  exponent  at  t = 15  second  is  the 
estimated  rx) . The  Lexp  can  be  calculated  by  the  following 
equations: 

M 

sum  = Z irexpi#  (3.4) 

i=l 

Lexp  = sum  / M,  (3.5) 

where  irexp-^  is  the  ith  instantaneous  running  exponent 
defined  in  equation  (3.1)  and  M is  the  number  of  iterates. 

In  this  study,  all  the  segments  are  categorized  as  the  stage 
two  sleep  EEG.  The  last  column  of  Table  3-2  lists  the  final 
estimate  of  T1  after  propagating  throughout  each  15-second 
segment.  To  describe  how  quickly  the  running  exponents 
converge.  Table  3-2  shows  the  average  and  standard  deviation 
of  the  running  exponents  for  the  last  half  and  one-fifth  of 
the  propagating  process.  Those  six  Lexp-t  curves  show 
convergence  to  within  a tolerable  range  of  deviation  (under 


(bits/  sec) 
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Figure  3-5  Running  exponent  vs.  time  curves  for  six 
15-second  sleep  stage  two  EEG  segments. 
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Table  3-2 

Fx  for  Six  Consecutive  15-Second  EEG  Segments 


segment 

50%a 

20%a 

ri 

l-15c 

1.57 

± 0.11 

1.47  ± 0.04 

1.43 

16-30 

1.83 

± 0.10 

1.84  ± 0.09 

2.00 

31-45 

1.76 

± 0.08 

1.69  ± 0.07 

1.55 

45-60 

1.91 

± 0.08 

1.86  ± 0.07 

1.89 

61-75 

1.79 

± 0.11 

1.86  ± 0.05 

1.83 

76-90 

2.31 

± 0.23 

2.09  ± 0.04 

2 . 10 

a)  n%:  the  percentage  of  the  total  running  exponents  used  to 
compute  the  average  and  standard  deviation,  starting 
from  the  last  running  exponent  (r^)  and  tracing 
backwards  to  involve  the  last  n%  of  the  entire 
running  exponents. 
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4.8%  of  the  average)  for  the  last  one-fifth  propagation. 
Unfortunately,  the  results  calculated  from  different 
segments  of  the  same  brain  state  are  spread  over  the  range 
1.43  to  2.10  bits/sec! 

The  inconsistency  of  the  estimated  largest  Lyapunov 
exponent  may  result  from  a change  in  the  system  parameters 
(see  explanation  in  Chapter  2 and  5). 

The  second  issue  which  may  cause  undesirable  effects  is 
the  noise-contamination  condition.  The  noise  component 
undoubtedly  is  detrimental  to  the  optimization  of  local 
structure  (16).  It  has  been  noted  that  Lyapunov  exponents 
are  not  well  defined  in  the  presence  of  any  amount  of  noise. 
In  (67,70),  it  caused  mostly  an  underestimate  of  Lyapunov 
exponents.  Moreover,  the  noise-contaminated  environment  may 
impact  negatively  the  determination  of  the  implementation 
parameters.  Therefore,  our  first  task  is  to  quantitatively 
characterize  the  noise  effect. 

Noise  Effect  on  T1  Estimation 

From  the  spectral  analysis  of  the  same  sleep  stage  II 
EEG,  two  peaks  located  at  49.80  and  70.31  Hz  are  identified 
to  be  noise  components,  and  the  peaks  at  2.44  and  14.16  Hz 
are  from  the  signal  itself.  First  of  all,  the  effect  of 
cutoff  frequency  on  the  estimated  r:  is  investigated.  The 
signal  (the  third  15-second  segment)  is  filtered  by  a 
Butterworth  lowpass  filter  (N  = 9)  with  cutoff  frequency  fc 
ranging  from  15  to  70  Hz.  We  previously  verified  that  the 
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correlation  dimension  is  affected  very  little  by  the  cutoff 
frequency  if  the  filter  does  not  distort  those  important 
patterns  of  a signal.  Accordingly,  it  is  not  necessary  to 
analyze  the  dimensionality  for  all  the  twelve  cases 
(different  fc)  . Instead,  we  only  evaluate  the  dimension  of 
the  signal  filtered  with  fc  = 60,  45,  and  30  Hz.  The 
sequence  of  computations  includes  (1)  filtering,  (2)  MI 
calculation,  (3)  dimension  estimate,  and  (4)  rx  estimate. 

The  MI  curves  are  plotted  in  Figure  3-6. 

Interestingly,  the  ambiguous  local  minimum  at  delay  = 2 
disappears  from  the  MI  curve  after  filtering  a considerable 
amount  of  noise.  All  the  minima  are  hardly  recognized. 

Hence  the  first  local  minimum  of  MI  curves  becomes  an 
appropriate  choice  for  delay.  From  the  analysis  of  the  CIM 
and  CEM  curves,  embedding  dimension  n = 6 is  selected  (/x  = 
4.625,  4.447,  and  4.268  for  fc  = 60,  45,  and  30  Hz, 
respectively) . With  all  the  other  parameters  unchanged 
(pscalmx  = 20%,  pscalmn  = 5%,  EVOLV  = 15),  the  estimation  of 
rl  and  the  values  obtained  by  averaging  the  final  one-fifth 
propagation  are  listed  in  Table  3-3  and  plotted  in  Figure  3- 
7.  The  time  delay  used  for  each  case  is  selected  from  the 
first  local  minimum  of  the  corresponding  MI  curve  (except 
the  case  fc  = 65  Hz,  the  second  local  minimum  is  used).  The 
dashed  line  (Figure  3-7)  depicts  the  estimated  result  of  the 
unfiltered  segment  (T1  = 1.55  bits/sec).  The  range,  average 
± standard  deviation,  for  the  final  one-fifth  propagation  is 


Mutual  Information  (bits) 


140 


delay  (samples) 


Figure  3-6  MI  vs.  delay  curves  for  a 15-second  EEG 

filtered  by  Butterworth  lowpass  filter  with 
fc  = 60  Hz  (stars  and  solid) , 

45  Hz  (stars  and  dashed) , and 
30  Hz  (circles  and  solid) . 
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Table  3-3 

Effect  of  Cutoff  Frequency  on  rx 


*c 

delay 

rx 

last  20% 

(Hz) 

(samples) 

(bits/s) 

(bits/s) 

15.00 

12 

2.09 

2.09  ± 0.04 

20.00 

10 

3.26 

3.31  ± 0.05 

25.00 

8 

3.09 

3.07  ±0.03 

30.00 

8 

3.61 

3.59  ± 0.05 

35.00 

11 

4.20 

4.21  ± 0.05 

40.00 

8 

3.22 

3.30  ± 0.06 

45.00 

9 

4.12 

4.08  ± 0.05 

50.00 

7 

3.03 

3.17  ± 0.05 

55.00 

8 

3.15 

3.11  ± 0.06 

60.00 

7 

3.24 

3.31  ± 0.06 

65.00 

10 

3.73 

3.69  ± 0.09 

70.00 

8 

3.46 

3.55  ± 0.12 

(bits/ sec) 
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Figure  3-7 


Effect  of  cutoff  frequency  on  the  largest 
Lyapunov  exponent. 
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described  by  a short  vertical  bar;  and  the  circle  represents 
the  estimated  rx  for  each  fc.  Obviously,  there  is  no 
monotonic  property  relating  the  cutoff  frequency  of  a 
lowpass  filter  to  the  estimated  rlf  yet  the  exponent 
significantly  increases  after  filtering  (about  twice  the 
exponent  of  the  unfiltered  segment) . An  undesirable  fact 
observed  is  the  large  deviation  between  the  estimated 
exponents  of  segments  with  considerably  similar  waveforms, 
e.g.,  the  exponents  of  the  segment  filtered  with  fc  = 25.0 
and  30.0  Hz  are  3.09  and  3.61  bits/ sec,  respectively. 
Therefore,  we  conclude  that  linear  filtering  biases  in  an 
inconsistent  manner  the  largest  Lyapunov  exponent  and  should 
be  used  with  caution. 

Prior  to  the  next  study  of  EVOLV,  we  shall  further  take 
a note  about  the  "convergent  capability"  of  rx  by 
investigating  the  largest  Lyapunov  exponent  of  the  EEG 
segment  (filtered  with  fc  = 30  Hz)  with  different  durations. 
Effect  of  EEG  Duration  on  rx  Estimation 

Table  3-4  illustrates  the  estimates  of  for  the  EEG 
segment  with  different  durations.  The  rx  for  a 15-second 
segment  (3.61  bits/sec)  is  very  close  to  the  one  for  a 60- 
second  segment  (3.22  bits/sec),  though  some  exponents 
estimated  "midway"  deviate  from  it  (such  as  Tj  = 2.82,  2.83, 
and  4.21  bits/sec  for  duration  of  35,  40,  and  20  seconds 
respectively) . Hereof  we  shall  recalculate  the  exponents 
for  the  six  15-second  segments  (Table  3-2)  after  filtered  by 
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Table  3-4 

Convergent  Capability  of  Estimated  rx 


for  Sleep 

Stage  Two  EEG 

duration 

(seconds) 

Delay 

(samples) 

Last  20% 
(bits/sec) 

rl 

(bits/sec) 

15. 0a 

8 

3.59  ± 0.05 

3.61 

20.0 

9 

4.30  ± 0.08 

4.21 

25.0 

10 

3.65  ± 0.04 

3.68 

30.0 

10 

3.68  ± 0.05 

3.63 

35.0 

15 

2.85  ± 0.03 

2.82 

40.0 

17 

2.88  ± 0.02 

2.83 

45.0 

17 

3.08  ± 0.03 

3.09 

50.0 

13 

3.02  ± 0.03 

3.09 

55.0 

13 

3.20  ± 0.08 

3.33 

60.0 

13 

3.29  ± 0.05 

3.22 

This  segment 

is  the  third 

one  of  the  six 

consecutive  15- 

second  segments,  as  the  waveform  in  this  portion  is  most 
representative  for  sleep  stage  two  EEG. 
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the  Butterworth  LPF,  fc  = 30  Hz,  since  the  results  may 
provide  the  information  for  explaining  the  varying  results 
in  Table  3-4.  Six  exponents  obtained  for  the  six 
consecutive  segments  are  4.61  (7),  3.70  (10), 3. 61  (8),  4.13 
(9),  3.01  (13),  and  5.14  (7)  respectively,  where  the  number 
in  parentheses  is  the  delay  in  samples.  The  first  four 
values  relate  somewhat  to  those  listed  in  Table  3-4.  The 
slowly  divergent  activities  during  the  3rd  EEG  segment  bias 
the  estimations  for  the  35-second  and  40-second  segment  and 
result  in  smaller  exponents  (2.82  and  2.83).  Then  the 
fiducial  trajectory  propagates  through  the  space  of  fast 
expansion  (1^  = 4.13  from  the  4th  segment)  and  consequently 
the  exponent  gradually  increases  (3.09  -+  3.09  -*•  3.33  -*• 

3.22).  The  above  experiment  illustrates  the  sensitive 
dependence  of  r2  on  the  signal  characteristics.  Due  to  the 
time  varying  nature  of  the  EEG,  we  might  foresee  that,  even 
with  a duration  » 60  seconds,  the  estimated  exponent  still 
varies  up  and  down  all  the  way  through  the  end  of  sleep  two 
stage. 

It  should  be  pointed  out  here  that,  while  searching  for 
the  way  of  relating  the  complicated  real  biological  system 
with  the  artificial  neural  network  through  the  analysis  of 
their  dynamical  properties,  we  need  to  work  it  out  under 
some  assumptions  (e.g.,  the  requirement  of  the  signal 
duration)  to  simplify  the  complexity  of  the  investigating 
circumstances.  At  this  point,  we  shall  not  waste  our  effort 
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on  seeking  a segment  length  for  a consistent  exponent 
estimate,  which  may  never  occurs  in  a dynamical  system 
having  time-varying  parameters.  Fortunately,  the  Lexp-t 
curve  (Figure  3-5)  for  a given  15-second  segment  is  able  to 
converge  satisfactorily  during  the  last  one-third 
propagation.  It  indicates  that  the  duration  is  sufficient 
to  allow  the  convergence  of  the  running  exponent.  The 
effect  of  insufficient  points  on  the  running  exponent  vs. 
time  curve  is  clearly  illustrated  in  Figure  3-8,  the 
temporal  convergence  of  running  exponent  for  segment  lengths 
of  4,  8,  15,  30,  and  60  seconds.  Obviously,  the  Lexp-t 
curves  for  the  4-second  and  8-second  segments  ( F-j^  = 4.30  and 
2.69  bits/sec)  still  oscillate  vigorously;  while  with  a 15- 
second  duration  the  running  exponent  settles  down  to  a 
tolerable  range.  Hence,  instead  of  struggling  throughout 
the  15000-point  attractor  which  takes  about  65  minutes  of 
computation  (pscalmx  = 20%,  pscalmn  = 5%,  EVOLV  = 15 
samples)  for  the  12  MHz  IBM-AT  with  coprocessor,  a segment 
length  of  15  seconds  (3750  samples)  will  be  used  to  study 
the  effect  of  other  parameters  in  the  estimation  of  the 
largest  Lyapunov  exponent  (the  computation  time  is  about  7 
minutes  with  these  parameters) . 

Evolution  Time  EVOLV  between  Replacement 

In  our  previous  study  of  the  signal  duration,  a fixed 
EVOLV  of  15  samples  was  chosen  according  to  the  following 
explanation.  From  the  power  spectrum,  the  signal's  energy 
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Figure  3-8  Running  exponent  vs.  time  curves  for  different  segment  lengths. 
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are  distributed  mostly  over  the  range  of  frequency  < 14  Hz, 
within  which  the  major  peak  of  spectrum  is  located  at  2.4414 
Hz.  With  a sampling  rate  of  250  Hz,  there  are  at  least  18 
samples  in  a period  of  a sinewave  of  14  Hz  or  less.  As 
mentioned  before,  the  "rhythm"  of  the  waveform  relates  to 
the  orbital  movement  of  the  reconstructed  attractor  and  the 
possibility  of  folding  or  stretching  process.  The  "18 
samples",  while  not  being  able  to  measure  the  orbital  size, 
can  provide  the  information  of  how  fast  the  skewness  occurs 
along  a certain  axis.  Based  on  the  suggestion  of  choosing 
an  EVOLV  in  the  range  0.5  to  1.5  orbits  (16),  we  selected 
the  EVOLV  of  15  samples  in  the  preliminary  study.  In  this 
section,  the  fixed  EVOLV  algorithm  is  implemented  to 
investigate  the  effect  of  EVOLV  on  the  exponent  estimate  and 
to  find  the  range  of  EVOLV  which  provides  a stable  estimate 
of  r^. 

The  sleep  stage  two  EEG  (N  = 3750)  is  first  filtered  by 
a Butterworth  lowpass  filter  (fc  = 30  Hz).  As  shown 
previously,  the  attractor  is  reconstructed  in  a six- 
dimensional space  with  time  delay  of  8 samples  (0.032 
second)  which  is  the  first  minimum  of  the  MI  curve.  The 
length  scale  of  the  local  structure  is  circumscribed  by 
pscalmx  = 20%  and  pscalmn  = 5%.  with  the  fixed  propagation 
step,  rx  is  evaluated  for  the  evolution  time  EVOLV  ranging 
from  l sample  (0.004  second)  to  41  samples  (0.164  second)  in 
steps  of  2 samples.  The  result  shown  in  Figure  3-9,  rx  vs. 
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EVOLV,  is  of  great  interest  since  it  displays  two  ranges  of 
EVOLV  that  provide  reasonably  stationary  exponent  estimates, 
rx  = 3.53  ± 0.20  (11  < EVOLV  < 21)  and  2.52  ± 0.14  (29  < 
EVOLV  < 39) . For  EVOLV  < 5,  large  exponents  are  obtained 
due  to  the  accumulation  of  orientation  errors  (16) . On  the 
other  hand,  the  exponent  significantly  drops  with  EVOLV  > 

23,  which  is  probably  caused  by  the  expanding  and 
contracting  activities  during  an  EVOLV-sample  propagation. 

Actually  we  have  found  that,  within  a typical  one-step 
evolutionary  local  structure  [P(t),Q(t')j  -*• 

[P (t+EVOLV) , Q (t 1 +EV0LV) ] , the  separation  vectors  (Figure  3- 
10)  between  the  fiducial  trajectory  P(t)  -*•  P (t+EVOLV)  and 
the  selected  nearby  trajectory  Q(t')  -*■  Q(t'+EVOLV)  (Figure 
3-1)  oscillate  in  magnitude  for  EVOLV  greater  than  20 
samples!  As  is  stressed  in  (16),  such  an  oscillation  either 
indicates  the  stretching  and  folding  process,  and  a 
replacement  for  Q is  needed;  or,  the  oscillation  in  distance 
(P/Q)  can  be  neglected  since  it  might  be  caused  by  a 
variation  of  the  attractor's  shape.  The  evolved  distance  of 
each  propagating  small-scaled  local  structure  (P,Q)  is 
displayed  in  Figure  3-10  for  the  distances  of  the  first  1000 
separation  vectors,  which  are  separated  into  small-scaled 
local  structures  (with  EVOLV  = ll,  21,  and  29)  by  the 
vertical  grids.  The  chart  exhibits  a vigorous  oscillation 
of  separation  vectors,  from  which  we  expect  the  dramatic 
complexity  of  the  attractor's  structure  and  possibly, 
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Figure  3-9  Estimated  largest  Lyapunov  exponent  for 
1 < EVOLV  < 41  samples. 


Distance  between  fiducial  and  qualified  nearby  trajectory 
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Figure  3-10  DPg  maps  for  filtered  sleep  stage  two  EEG  (f  =30  Hz)  . 

(a)  EVOLV  = 11;  (b)  EVOLV  = 21;  (c)  EVOLV  = 29. 


Distance  between  fiducial  and  qualified  nearby  trajectory 
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frequent  folding  and  stretching  activities.  Consequently, 
these  facts  support  the  choice  of  the  range  11  < EVOLV  < 21 
as  appropriate  for  the  computation  of  the  largest  Lyapunov 
exponent . 

The  horizontal  dashed  line  in  Figure  3-10  and  3-11 
represents  the  maximum  and  minimum  scale  limitations  (SCALMX 
and  SCALMN)  of  each  local  structure.  A separation  beyond 
SCALMX  indicates  a replacement  for  Q is  needed.  The 
magnitude  of  a separation  vector  will  be  called  DPQ 
hereinafter  as  it  is  equivalent  to  the  distance  between  a 
point  P(t+ts)  on  the  fiducial  trajectory  and  another  point 
Q(t'+ts)  on  the  selected  nearby  trajectory  (ts  < EVOLV).  In 
the  following  we  use  the  Lorenz  model  (Figure  3-4)  to 
describe  the  usefulness  of  the  DPq  map  to  help  us  visualize 
the  effect  of  these  parameters  in  the  exponent  estimation. 
Figure  3-11 (a)  shows  an  example  of  DPq  map  for  an  EVOLV  that 
is  too  small.  The  fluctuation  is  about  half  a cycle  within 
a local  structure  defined  by  EVOLV  = 9,  which  leads  to  an 
overestimate  of  r1  (2.46).  On  the  other  hand,  the 
oscillations  of  DPQ  for  EVOLV  = 19  (Figure  3-11 (b) ) or 
larger  show  the  folding  and  stretching  effects  (large 
oscillation),  which  will  underestimate  the  exponent  (1.79). 
With  EVOLV  < 13  samples,  is  slightly  overestimated  as 
shown  in  the  following: 

EVOLV 5 7 9 11 13 15 17 19 

ri  2.07  2.28  2.46  2.46  2.53  2.16  2.11  1.79 
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Figure  3-11  DPg  maps  for  Lorenz  attractor  (pscalmx=30% , pscalmn=5%) . 
(a)  EVOLV  =9;  (b)  EVOLV  = 19. 
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Comparing  Figure  3-10  and  3-11,  we  notice  the  DpQ  for  EEG 
varies  tremendously  in  magnitude,  which  indicates  that  there 
are  many  small-scale  transitions  on  the  EEG  attractor's 
structure.  We  wonder  if  we  should  closely  trace  each 
transition  (i.e.,  a cycle)  in  a DpQ  map.  If  yes,  according 
to  Figure  3-10,  we  should  not  use  an  EVOLV  larger  than  20 
samples  for  the  sleep  stage  two  EEG  with  a sampling  rate  of 
250  Hz.  Otherwise,  a large  orientation  error  may  result 
from  a small  EVOLV,  which  causes  an  overestimate  of  1^. 
However,  it  still  remains  guestionable  whether  the 
oscillation  of  separation  vectors  results  from  a stretching 
and  folding,  or  is  simply  an  indication  of  a variation  of 
the  attractor's  shape.  Since  the  range  9 < EVOLV  < 19 
(Table  3-5)  provides  stable  exponent  estimates  for  sleep 
stage  two  EEG,  we  will  use  an  EVOLV  of  15  samples  for  the 
study  of  SCALMX  and  SCALMN. 

Length  Scales  Limitation  (SCALMX  and  SCALMN) 

The  use  of  pscalmx  =20%  and  pscalmn  = 5%  is  based  on 
some  published  works  (8,16).  However  these  parameters  are 
highly  data-dependent . The  DPq  map  introduced  in  previous 
section  provides  some  useful  information.  Later  we  use  a 
"frame"  to  represent  the  magnitude  of  separation  vectors  in 
each  EVOLV-sample  span  fenced  by  the  gridding  lines  as  shown 
in  Figure  3-10  to  3-12. 

°f  all,  it  is  obvious  that  the  "instantaneous 
running  exponent"  (irexp  defined  in  equation  (3.1))  can  be 
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Table  3-5 

rx  Estimated  for  Different  EVOLV,  pscalmx,  and  pscalmn 


(pscalmx=20%)  (EV0LV=15)  (pscalmx=20%) 

(pscalmn=5%)  (pscalmn=5%)  (EV0LV=15) 

EVOLV  r2  pscalmx  r2  pscalmn  rx 

(samples)  (bits/s)  (%)  (bits/s)  (%)  (bits/s) 


1 

7.38 

10 

3.74 

1 

2.67 

3 

5.19 

12 

4.49 

2 

2.67 

5 

4.97 

14 

3.71 

3 

2.67 

7 

4.01 

16 

3.46 

4 

3.68 

9 

3.87 

18 

4.28 

5 

3 . 61 

11 

3.80 

20 

3.61 

6 

3.59 

13 

3.50 

22 

2.96 

7 

3.41 

15 

3.61 

24 

2.24 

8 

3 .40 

17 

3.28 

26 

1.99 

9 

3.40 

19 

3.27 

28 

1.62 

10 

2.95 

21 

3.71 

30 

1.63 

11 

3 . 56 

23 

2.81 

32 

1.04 

12 

3 . 39 

25 

2.43 

34 

0.92 

13 

2 .38 

27 

3.04 

36 

0.67 

14 

2 . 73 

29 

2.41 

38 

0.91 

15 

2 . 59 

31 

2.65 

40 

0.72 

33  2.37 


35  2.61 
37  2.37 
39  2.72 
41  2.18 
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determined  from  the  DPq  map.  Since  the  instantaneous 
running  exponent  is  proportional  to  the  logarithm  of  the 
fraction  df  / di  (di  and  df  are  the  initial  and  final 
separation  of  the  segmental  fiducial  trajectory  and  a 
selected  nearby  trajectory  segment) , the  beginning  and 
ending  value  of  a DPQ  segment  within  a certain  frame  (Figure 
3-10  and  3-11)  determines  the  instantaneous  running  exponent 
at  that  moment.  Moreover,  the  DpQ  map  reveals  how  the 
expanding  and  contracting  activities  (i.e.,  the  positive  and 
negative  instantaneous  running  exponents)  are  distributed 
over  the  entire  attractor.  For  instance,  the  DpQ  maps  of 
Lorenz  attractor  (Figure  3-11  and  3-12)  show  an  abrupt 
orbital  divergence  after  the  trajectory  evolving  for  about  5 
seconds  (500  samples) . it  indicates  that  the  local 
structure  chosen  at  that  time  is  very  likely  to  involve 
points  which  fly  to  opposite  lobes.  The  result  is  an 
enormous  contribution  to  the  exponent  estimate  (16) . 

Secondly,  from  the  Dpg  map,  we  know  when  a replacement 
is  taking  place  and  how  many  times  replacement  occurs  during 
the  estimating  process.  At  the  end  of  each  frame,  if  the 
Dpq  goes  above  the  dashed  line  (SCALMX) , a replacement  is 
needed  and  there  exists  the  discontinuity  from  here  to  the 
beginning  of  next  frame.  However,  if  the  DPQ  overshoots 
SCALMX  at  any  other  time,  there  will  be  no  inspection  and 
accordingly,  no  replacement.  The  DpQ  map  gives  a panoramic 
view  to  help  us  understand  the  role  of  EVOLV  and  SCALMX  in 
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the  process  of  estimating  the  largest  Lyapunov  exponent.  In 
addition,  from  the  DpQ  map  we  may  have  a more  concrete  idea 
of  how  the  EVOLV  and  SCALMX  work  together  to  determine  the 
formation  of  each  evolutionary  local  structure. 

Observing  the  10-second  DpQ  maps  of  the  Lorenz 
attractor  in  Figure  3-12  where  the  parameters  are  all  the 
same  (pscalmn  = 5%  and  EVOLV  = 15  samples)  with  the 
exception  of  pscalmx  = 30%  (Figure  3-12 (a))  and  35%  (Figure 
3-12 (b) ) , we  may  get  a picture  of  how  the  SCALMX  affects  the 
exponent  estimate.  The  running  exponents  of  both  cases  are 
the  same  until  propagating  for  3.15  seconds  (315  samples); 
then  at  the  instant  marked  by  a cross,  the  magnitude  of  DPQ 
comes  between  30%  and  35%  of  the  spatial  extent,  which 
results  in  two  completely  different  DPQ  maps  (in  Figure  3- 
12 (a)  a replacement  occurs,  while  in  Figure  3-12 (b)  no 
replacement  occurs).  The  estimated  r1  is  2.16  (pscalmx  = 

30%)  and  1.76  (pscalmx  = 35%).  Apparently,  some  small 
variation  in  the  implementing  parameter  may  lead  to  very 
different  results.  In  addition,  we  have  noted  that  the  same 
xponent  can  be  obtained  with  an  appropriate  coordination  of 
EVOLV  and  SCALMX  (pscalmx) . Taking  the  Lorenz  attractor 
(Figure  3-4)  as  an  example,  we  tried  to  vary  both  the  EVOLV 
and  pscalmx  (pscalmn  = 5%)  to  get  an  accurate  estimate. 
Interestingly,  the  estimate  is  close  to  the  one  reported  in 
( ) (2.16  bits/ sec)  with  several  combinations:  (1)  EVOLV  = 

16,  pscalmx  = 28%  (2.15),  (2)  EVOLV  = 15,  pscalmx  = 30% 
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(2.16),  (3)  EVOLV  = 14,  pscalmx  = 32%  (2.17),  and  (4)  EVOLV 

= 13,  pscalmx  = 35%  (2.18). 

From  the  above  experiment  and  the  results  in  Table  3-5, 
we  find  that  increasing  the  SCALMX  has  the  same  effect  as 
increasing  the  EVOLV.  Both  result  in  a smaller  exponent. 

Our  study  of  the  sleep  stage  II  EEG  shows  that  approximately 
stationary  exponent  estimates  can  be  obtained  with  10%  < 
pscalmx  <20%  when  EVOLV  = 15  samples. 

Since  short  separation  vectors  are  likely  to  involve 
large  amount  of  noise,  the  SCALMN  (pscalmn)  is  employed  to 
counteract  this  problem,  which  may  be  still  significant 
after  filtering.  Unfortunately,  we  could  not  see  a clear 
way  to  derive  the  value  of  SCALMN  from  the  DPQ  map,  the 
signal  itself,  or  the  phase  trajectory.  Nevertheless  our 
previous  experiences  in  the  evaluation  of  EEG's  correlation 
dimension  inspire  us  to  think  of  using  the  "knee"  phenomenon 
as  a criterion.  The  CIM  curves  for  the  segment  is  shown  in 
Figure  3-13  for  embedding  dimension  from  2 to  18. 

Inspecting  the  CIM  curve  for  n = 6,  a circular  mark  at  ln(r) 
— 4.0  divides  the  curve  into  two  portions:  the  upper  region 
characterizing  the  deterministic  EEG  and  the  lower  region 
where  the  noise  activities  settle  in.  The  logarithm  of  the 
spatial  extent  of  the  EEG  segment  is  7.0424,  in  other  words, 
the  largest  separation  between  any  pair  of  the  points  in  the 
embedding  space  is  1144.  From  the  CIM  curves  we  see  that 
the  intermediate  knee  occurs  when  r = In'1 (4.0)  = 54.5. 


ln(C(r)) 


161 


Figure  3-13  CIM  curves  for  filtered  sleep  stage  two  EEG. 
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Therefore,  the  noise  scale  starts  at  4.77%  of  the  spatial 
extent.  In  the  right  column  of  Table  3-5,  the  exponents  are 
stable  in  two  ranges  of  pscalmn,  1%  < pscalmn  < 3%  and  4%  < 
pscalmn  < 9%.  However  for  pscalmx  < 3%,  noise  contamination 
still  partially  affects  the  estimation  and  causes  an 
underestimate  (2.67  bits/sec)  as  which  observed  in  the 
exponent  estimates  for  the  unfiltered  segments  (Figure  3-5 
and  Table  3-2) . 

Lyapunov  Exponent  of  a Spike  and  Wave  Seizure 

It  is  worth  estimating  the  largest  Lyapunov  exponent  of 
a more  regular  signal  pattern  (such  as  the  EEG  in  a spike 
and  wave  discharge)  since  the  sleep  stage  two  EEG  involves 
more  degrees  of  freedom  and  lots  of  mixed  waves.  These 
factors  complicate  the  exponent  estimation  and  determination 
of  implementation  parameters. 

The  procedure  and  result  of  estimating  r2  for  the  EEG 
with  epileptic  spikes  (Figure  3-14)  are  presented  below  to 
conclude  this  chapter.  First  of  all,  the  CIM  curves  (Figure 
3-15)  and  MI  curves  provide  information  of  selecting  n and 
r.  They  are  n = 5 and  r = 8 respectively.  From  the  power 
spectrum  (Figure  3-16)  50%  of  the  total  power  is  distributed 
over  the  frequency  range  0.0  Hz  < frequency  < 6.59  Hz.  The 
EVOLV  will  be  determined  according  to  the  frequency  (fH  = 

6.59  Hz)  that  equally  divides  the  total  power.  Hence  the 
EVOLV  around  38  samples  is  investigated.  The  minimum  scale 
limitation  (SCALMN)  can  be  computed  from  the  correlation 
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Figure  3-14  EEG  segment  with  epileptic 
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Figure  3-16  Power  spectrum  for  EEG  segment  with 
epileptic  spikes. 
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integral  curve  (n  = 5) . On  the  curve  a circular  mark  at 

ln(r)  = 5.1  approximately  separates  the  regions 

characterizing  the  signal  and  noise  scale.  So  SCALMN  « ln“ 

1 ( 5 . 1)  = 164.022  is  used.  Spatial  extent  (see  definition  in 

Appendix)  can  also  be  roughly  estimated  from  the  CIM  curves. 

For  example,  ln[C(r)]  = 0.0  for  ln(r)  > 8.3  (n  = 5)  and 

ln[C(r) ] + 0.0  for  ln(r)  < 8.1  (n  = 5),  therefore  3294.47  < 

extent  < 4023.87  (extent  = 3604.58).  Then  pscalmn  « 4.55% 

is  obtained.  It  takes  some  effort  to  determine  the  last 

parameter  (SCALMX  and  pscalmx) . No  concrete  criterion  is 

found  for  selecting  the  maximum  scale  limitation  of  EEG.  At 

♦ 

this  stage,  we  propose  to  check  for  exponent  stability  over 
a wide  range  of  pscalmx.  Figure  3-17  plots  the  estimated 
largest  Lyapunov  exponent  for  pscalmx  from  10%  to  40%.  It 
is  evident  that  the  estimation  displays  a plateau  for  10%  < 
pscalmx  < 25%.  Then  the  final  estimation  of  (4.20  ± 

0.34)  is  obtained  by  averaging  the  exponents  within  this 
pscalmx  range  (discarding  the  two  largest  and  smallest 
exponents  when  pscalmx  = 11%,  12%,  18%,  and  21%) . The 
dashed  horizontal  line  in  Figure  3-17  shows  the  average  of 
the  stable  exponents.  This  approach  cannot  be  applied  to 
the  selection  of  EVOLV  and  pscalmn  since  more  than  one  range 
of  stationary  estimation  are  found  for  the  r2  vs.  EVOLV  and 
rl  vs.  pscalmn  (Figure  3-9  and  Table  3-5) . 
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Figure  3-17  Largest  Lyapunov  exponent  vs.  pscalmx 

(N=3000 , n=5 , 7=8,  EVOLV=38,  pscalmn=4 . 55%) . 
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Discussion  and  Conclusion 

We  have  shown  our  systematic  approach  of  quantifying 
the  effect  of  the  implementing  parameters  used  in  the 
estimation  of  the  largest  Lyapunov  exponent  for  experimental 
data.  A positive  rx  (3.61  bits/sec)  is  estimated  for  sleep 
stage  two  EEG  segment  with  N = 3750,  n = 6,  delay  = 8,  EVOLV 
= 15,  pscalmx  = 20%,  and  pscalmn  = 5%.  This  result  further 
corroborates  the  correlation  dimension  evidence  that  sleep 
stage  two  EEG  segment  possesses  an  attractor  that  is 
chaotic.  However  the  actual  value  of  the  largest  Lyapunov 
exponent  can  not  be  resolved.  Our  estimate  is  quite 
different  from  the  result  reported  in  (6)  (1^  between  0.4 

and  0.8  for  sleep  stage  two).  As  demonstrated  in  the  above 
experiments,  due  to  the  extreme  sensitivity  of 
estimation,  even  a slight  difference  in  the  implementing 
parameters  can  cause  a great  deviation  in  estimates.  In 
addition,  brain  function  is  assumed  time-varying.  Under  the 
circumstances  there  exist  even  more  difficulties  in  reaching 
a close  estimate  especially  among  laboratories  having 
different  sources  of  signals. 

In  this  Chapter,  we  have  proposed  our  procedure  of 
estimating  the  largest  Lyapunov  exponent  for  EEG  signals. 

In  summary,  the  EVOLV  is  determined  first  from  the  following 
equation: 


\ 


EVOLV  = fs  / fp, 


(3.6) 
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where  fg  is  the  sampling  frequency  and  fp  is  the  frequency 
component  that  equally  divides  the  total  power  of  the 
signal.  Then  extent  (the  largest  interpoint  distance  in  the 
trajectory)  is  computed.  To  select  the  minimum  scale 
limitation  (SCALMN  and  pscalmn) , we  propose  to  use  the 
information  of  CIM  curves  ("knee"  phenomenon) . The 
remaining  parameter,  SCALMX,  is  determined  as  follows.  The 
largest  Lyapunov  exponent  is  estimated  for  a wide  range  of 
pscalmx  so  that  the  range  of  stationary  estimation  can  be 
found.  The  final  result  is  then  computed  by  averaging  the 
exponents  within  the  range  providing  stationary  estimation. 

From  this  study,  we  have  made  a few  constructive 
propositions  for  the  exponent  estimation.  First  of  all,  we 
introduce  the  concept  of  using  information  of  a signal  (such 
as  power  spectrum)  to  estimate  the  possible  "orbital  size" 
and  the  appropriate  EVOLV.  And  we  propose  to  use  the  knee 
phenomenon  in  CIM  curves  to  estimate  the  minimum  length 
scale  limitation  (SCALMN) . We  have  also  shown  that  the  DPQ 
map  potentially  provides  useful  information  for  further 
investigation  of  local  structure  activity  of  EEG 
trajectories.  This  helps  us  develop  a structure-dependent 
algorithm  with  variable  EVOLV  to  estimate  the  largest 
Lyapunov  exponent  for  the  complex  EEG  trajectory. 


CHAPTER  4 

DIMENSION  AND  EXPONENT  ESTIMATES 
FOR  OTHER  BRAIN-STATE  EEG  SEGMENTS 

The  origin  of  our  study  is  based  on  the  hypothesis  that 
the  chaotic  nature  of  EEG  attractors  results  in  the 
capability  of  the  brain  to  generate  and  process  information. 
Concretely  speaking,  the  estimation  of  correlation  dimension 
allows  us  to  characterize  the  degree  of  complexity  of  the 
ongoing  EEG  dynamics  (an  increase  in  the  dimension 
corresponds  to  an  increase  in  mental  activity) . The 
evaluation  of  Lyapunov  exponent  provides  a possibility  of 
characterizing  the  local  structure  of  EEG  attractors.  Our 
methodology  developed  previously  has  been  applied  to  the 
sleep  stage  two  EEG  and  shown  that  it  possesses  the 
characteristics  of  low-dimensionality  and  chaotic  behavior 
(the  estimated  largest  Lyapunov  exponent  is  positive) . At 
this  point,  it  is  of  paramount  importance  to  explore  and 
compare  the  dynamical  properties  of  other  brain  states  to 
prove  or  disprove  our  hypothesis. 

In  this  study,  we  will  display  the  results  of  analyzing 
the  segments  under  sleep  stage  zero,  two,  three,  five 
(expressed  as  0,  II,  III,  V later),  EEG  with  epileptic 
spikes,  EEG  contains  muscle  activity,  and  awake  EEG  (eyes 
open)  with  spike-and-wave.  The  epileptic-spike  EEG  and 
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awake  EEG  with  spike-and-wave  are  from  one  subject;  the 
others  are  recorded  from  another  subject.  The  main  skeleton 
of  this  section  is  to  quantify  the  topological  properties  of 
the  EEG  attractors  by  means  of  dimensional  analysis  and 
exponent  estimation. 

Correlation  Dimension 

The  estimation  of  correlation  dimension  involves  two 
steps.  First  of  all,  the  minimum  number  of  data  points 
required  for  obtaining  a convergent  estimate  must  be 
determined  for  different  brain-state  EEG  segments.  At  this 
stage  N = 2000  (only  for  the  epileptic-spike  segments) , 

3000,  4000,  5000,  and  6000  samples  are  used.  For  noise- 
contaminated  signals,  the  original  and  filtered  segments  are 
analyzed  to  quantify  the  noise  effect.  For  each  N-sample 
segment,  the  mutual  information  needs  to  be  recalculated  to 
find  the  appropriate  time  delay.  Then  the  ln[C(r)]  vs. 
ln(r)  will  be  calculated  for  the  embedding  dimension  from  2 
to  20.  After  determining  all  the  implementation  parameters 
for  each  brain-state  EEG,  correlation  dimension  will  be 
calculated  for  four  different  segments  within  the  same  brain 
state. 

Figure  4-1  to  4-7  display  the  raw  EEG  segments  for  the 
seven  brain  states  (sampling  rate  = 250  Hz,  8-bit  ADC,  4 
channels) . The  sleep  EEGs  used  in  this  study  are  recorded 
from  CZ-A2 , and  EEG  segments  with  epileptic  spikes  are  from 
F3-A1.  The  segments  of  sleep  stage  0,  III,  and  V are  noisy. 


SPIKE.DAT  (channel  1,  1-3200  samples) 
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Figure  4-1  EEG  segment  with  epileptic  spikes. 


SLEEP1.DAT  (channel  1,  1-4000  samples) 
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Figure  4-2  Sleep  stage  0 EEG. 


SLEEP  STAGE  II  EEG  (channel  3,  i— 4000  samples) 
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Figure  4-3  Sleep  stage  II  EEG. 


SLEEP3.DAT  (channel  1,  1—4000  samples) 
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Figure  4-4  Sleep  stage  III  EEG. 


SLEEP5.DAT  (channel  1,  1—4000  samples) 
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Figure  4-5  Sleep  stage  V EEG . 


EEG  contains  Muscle  Activity  (channel  3,  1-4000  samples) 
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Figure  4-6  EEG  contains  muscle  activity. 


Awake  EEG  with  Spikes-and— Wave  (channel  3,  1-4000  samples) 
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Figure  4-7  Awake  EEG  with  spikes-and-wave 
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A Butterworth  lowpass  filter  with  cutoff  frequency  = 45  Hz 
will  be  applied.  The  implementation  parameters  and 
estimated  correlation  dimensions  are  listed  in  Table  4-1. 

The  fourth  column  (n's)  lists  the  numbers  of  embedding 
dimension  from  which  the  correlation  dimension  is  derived  by 
averaging  the  corresponding  correlation  exponents.  These 
values,  obtained  when  parallelism  develops  in  the  CIM  curves 
until  the  knee  effect  settles  in  (68) , signify  the  possible 
range  of  embedding  dimension  of  the  reconstructed 
attractors.  Figure  4-8  to  4-14  plot  the  CEM  curves  for  the 
signals  shown  in  Figure  4-1  to  4-7.  In  Figure  4-8,  4-9,  4- 
11,  and  4-12,  the  CEM  curves  for  raw  (filtered)  signals  are 
plotted  as  solid  (dashed)  curves,  and  the  computed  data  are 
marked  by  the  centered  symbols  as  diamonds  (N  = 2000) , 
circles  (N  = 3000) , stars  (N  = 4000) , and  crosses  (N  = 

5000) . The  CEM  curve  for  N = 6000  nearly  coincides  with  the 
one  for  N = 5000,  and  therefore  is  not  plotted.  Several 
features  described  below  are  observed  from  the  results  shown 
in  Table  4-1,  the  CEM  curves,  and  the  information  of 
signals. 

(1)  Signals  having  consistently  sluggish  activities 
usually  possess  the  attractor  of  lower  dimensionality.  Our 
investigation  makes  us  infer  a possible  connection  between 
the  signal  and  dimensional  estimation — a signal  having  its 
power  dispersed  through  a wide  range  of  frequency  tends  to 
exhibit  higher  dimensionality.  For  instance,  correlation 


Table  4-1 

Dimension  Estimates  for  Different  EEG  Segments 


signal 

N delay  n's 

correlation 

dimension 

error 

SPIKE (R) a 

2000 

8 

4.5.6 

3.338 

+ 

0.128 

8.40% 

SPIKE (R} 

3000 

8 

4.5.6 

3.544 

+ 

0.173 

2.77% 

SPIKE (R} 

4000 

8 

4.5.6 

3.660 

+ 

0.215 

-0.42% 

SPIKE (R} 

5000 

8 

4,5.6 

3.729 

+ 

0.218 

-2.31% 

SPIKE m 

6000 

8 

4,5.6 

3 . 645 

+ 

0.256 

0.00% 

SPIKEmb 

2000 

9 

4,5,6 

3.305 

+ 

0.264 

13.77% 

SPIKE  CF} 

3000 

9 

4,5.6 

3.563 

+ 

0.288 

7.04% 

SPIKE  fF} 

4000 

9 

4,5.6 

3.718 

+ 

0.334 

3.01% 

SPIKE (F} 

5000 

9 

4.5.6 

3.800 

+ 

0.362 

0.86% 

SPIKE (F} 

6000 

10 

4.5,6 

3.833 

0.409 

0.00% 

SLEEPO  (R) 

3000 

5 

10.11.12 

7.355 

+ 

0.081 

9.35% 

SLEEPO fR} 

4000 

7 

10.11.12 

7.924 

+ 

0.154 

2.33% 

SLEEPO fRi 

5000 

7 

10.11.12 

8.113 

+ 

0.251 

0.00% 

SLEEPO  CF} 

3000 

5 

6,7 

4.697 

+ 

0.106 

7.63% 

SLEEPO  f F^ 

4000 

6 

6,7 

4.828 

+ 

0.099 

5.06% 

SLEEPO  CF} 

5000 

5 

6.7 

5.036 

+ 

0.142 

0.96% 

SLEEPO (F} 

6000 

5 

6,7 

5.085 

+ 

0.159 

0.00% 

SLEEPII fRl 

3000 

9 

6,7.8 

4.518 

+ 

0.201 

8.14% 

SLEEPII fR} 

4000 

9 

6.7.8 

4.770 

+ 

0.252 

3 . 01% 

SLEEPII (R1 

5000 

9 

6,7.8 

4.843 

+ 

0.266 

1.52% 

SLEEPII (R> 

6000 

9 

6,7.8 

4.918 

+ 

0.269 

0.00% 

SLEEPII fR} 

3000 

9 

—6,7,8 

4.523 

+ 

0.240 

— 

SLEEPIIim 

3000 

10 

7,8,9 

5.828 

+ 

0.082 

3.23% 

SLEEPIIIfR} 

4000 

11 

7.8.9 

5.916 

+ 

0.036 

1.76% 

SLEEPIII fR} 

5000 

11 

7,8.9 

6.022 

+ 

0.080 

0.00% 

SLEEPIIim 

3000 

10 

4,5.6 

3.889 

+ 

0.214 

7.22% 

SLEEPIII m 

4000 

10 

4,5.6 

3.982 

+ 

0.256 

5.00% 

SLEEPIII (F^ 

5000 

10 

4.5.6 

4.110 

+ 

0.276 

1.96% 

SLEEPIII (F^ 

6000 

10 

4,5,6 

4.192 

+ 

0.313 

0.00% 

SLEEPVm 

3000 

8 

11.12. 13 

8.650 

+ 

0.183 

2.28% 

SLEEPV(R) 

4000 

9 

11.12.13 

8.825 

+ 

0.322 

0.29% 

SLEEPVrRl 

5000 

9 

11.12.13 

8.851 

+ 

0.345 

0.00% 

SLEEPVm 

3000 

8 

5,6 

4.271 

+ 

0.067 

5.13% 

SLEEPVm 

4000 

8 

5.6 

4.375 

+ 

0.113 

2 .83% 

SLEEPVm 

5000 

8 

5.6 

4.474 

+ 

0.141 

0.62% 

SLEEPY m 

6000 

8 

5.6 

4.502 

+ 

0.137 

0.00% 
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Table  4-1 
(continued) 


signal 

N 

delay  n*s 

correlation 

dimension 

error 

MUSCLE (R} 

3000 

6 

10.11.12 

8.947 

+ 

0.572 

-23.70% 

MUSCLE (R) 

4000 

6 

10.11.12 

6.952 

+ 

0.161 

3.88% 

MUSCLE (R^ 

5000 

6 

10.11.12 

7.125 

+ 

0.175 

1.50% 

MUSCLE (R) 

6000 

6 

10.11.12 

7.233 

+ 

0.368 

0.00% 

AWSPWcm 

3000 

17 

8.9 

5.048 

+ 

0.000 

10.88% 

AWSPW (R\ 

4000 

17 

8.9 

5.459 

+ 

0.032 

3.61% 

AWSPW (R} 

5000 

17 

8,9 

5.670 

+ 

0.042 

-0.10% 

AWSPW m 

6000 

19 

8.9 

5.664 

+ 

0.035 

0.00% 

a)  R:  raw  EEG  segments. 

b)  F:  EEG  signals  filtered  by  Butterworth  LPF  (fc  = 45  Hz). 

c)  AWSPW:  awake  EEG  with  spike-and  wave. 


correlat ion  exponent 
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Figure  4-8 


CEM  curves  for  EEG  segment  with  epileptic 
spikes. 


correlation  exponent 
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Figure  4-9 


CEM  curves  for  sleep  stage  0 EEG 


correlation  exponent 
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Figure  4-10  CEM  curves  for  sleep  stage  TI  EEG 


correlation  exponent 
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Figure  4-11  CEM  curves  for  sleep  stage  III  EEG 


correlation  exponent 
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Figure  4-12 


CEM  curves  for  sleep  stage  V EEG. 


correlation  exponent 
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Figure  4-13 


CEM  curves  for  EEG  contains  muscle  activity. 


correlation  exponent 
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embedding  dimension  (n) 


Figure  4-14 


CEM  curves  for  awake  EEG  with  spikes-and-wake. 
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dimensions  of  (filtered)  sleep  EEG  show  that  /x(0)  > n(V)  > 
M(HI) ; and  from  the  spectrum  curves,  the  peak  of  power 
spectrum  of  sleep  stage  0 (III,  V)  EEG  is  located  at  9.2773 
(2.9297,  7.8125)  Hz,  and  11.13%  (5.27%,  10.74%)  of  the 
frequency  components  (frequency  resolution  = 0.24414  Hz)  has 
the  power  greater  than  one  tenth  of  the  maximum  power. 
Another  information  obtained  from  this  study  is  the 
increasing  coherence  in  EEG  dynamics  when  sleep  stage 
proceeds  through  the  deep  sleep  (6) . 

(2)  A smaller  delay  is  needed  for  a vigorously 
oscillating  waveform  (e.g.  sleep  stage  0 EEG  and  EEG 
contains  muscle  activity) . Generally  the  first  local 
minimum  of  the  MI  curve  shifts  towards  a smaller  value  of 
delay  after  the  reduction  of  noise  components  by  filtering. 

(3)  The  delay  used  for  phase  space  reconstruction  tends 
to  be  larger  for  the  EEG  attractor  having  lower 
dimensionality. 

(4)  The  last  column  in  Table  4-1  lists  the  deviation  in 
% from  the  desired  result  estimated  from  the  longest  data 
sequence  (N  = 5000  and  6000  for  an  original  sequence  and  the 
filtered  version,  respectively) . Apparently,  the 
correlation  dimension  for  EEG  segments  with  spikes  quickly 
converges  (deviation  < 3%)  with  a short  duration  of  N = 

3000.  For  the  other  EEG's,  the  estimated  correlation 
dimension  shows  small  deviation  (less  than  5.1%)  when  N > 
4000.  Even  for  the  high  dimensional  muscle  activity,  the 
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segment  length  of  4000  samples  already  provides  a convergent 
result.  Therefore,  it  is  appropriate  to  use  a segment 
length  of  3000  samples  for  the  epileptic-spike  EEG  and  4000 
samples  for  the  others  analyzed  in  this  study. 

(5)  For  the  epileptic-spike  EEG,  filtering  does  not 
decrease  the  correlation  dimension  (Figure  4-8)  as  presented 
in  Chapter  2 . This  is  because  the  percentage  of  power  above 
45  Hz  (cutoff  frequency  of  the  filter)  is  less  than  0.5%.  A 
further  investigation  shows  that  the  filtering  effect  on 
correlation  dimension  is  negligible  for  the  signal  having 
insignificant  amount  of  power  above  the  cutoff  frequency  of 
the  lowpass  filter.  On  the  other  hand,  a large  reduction  of 
the  dimension  is  observed  for  the  signal  having  more  power 
above  the  cutoff  frequency.  The  power  spectrum  of  sleep 
stage  V EEG  shows  that  more  than  14.5%  power  is  distributed 
in  the  frequency  range  f > 45  Hz.  And  a reduction  of  4.4 
for  the  correlation  dimension  and  6 for  the  embedding 
dimension  (n)  is  found  after  filtering.  Another  interesting 
phenomenon  for  this  EEG  segment  (Figure  4-5)  without 
filtering  is  the  disappearance  of  "knees"  in  CIM  curves 
(Figure  4-15) . 

(6)  The  EEG  segment  containing  muscle  activity  exhibits 
the  highest  embedding  and  correlation  dimension  (n  > 10,  /i  « 
7).  A big  deviation  is  observed  for  N = 3000  samples.  The 
29%  increase  in  correlation  dimension  (compared  to  the 
result  for  N = 6000)  is  probably  due  to  the  volume  expansion 
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Figure  4-15  CIM  curves  for  sleep  stage  V EEG. 
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of  the  trajectory  in  the  interval  between  the  2001st  and  the 
3000th  sample  (Figure  4-6,  the  third  trace  from  the  top) . 
Afterwards  the  EEG  activity  returns  to  its  usual  pattern. 

The  volume  of  the  trajectory  is  somewhat  reduced  which 
results  in  a smaller  correlation  dimension  when  4000  samples 
are  involved  in  the  estimation.  Here  we  conclude  that  the 
quantified  dynamical  characteristics  (e.g.  correlation 
dimension)  may  help  us  detect  the  change  of  activity  of  the 
original  one-dimensional  time  series  for  the  same  subject. 

(7)  It  should  be  noted  that  CEM  curves  for  the  filtered 
sleep  stage  0 and  III  EEG  segments  are  found  to  have  the 
tendency  to  merge  into  CEM  curves  for  the  raw  signals  when 
embedding  dimension  goes  higher  (in  Figure  4-6  and  4-7,  n > 
13  for  sleep  stage  0 and  n > 9 for  sleep  stage  III) . The 
only  explanation  we  may  conjecture  at  this  moment  is  that 
increasing  the  embedding  dimension  far  above  what  is  needed 
will  increase  the  level  of  noise-contamination  and  distort 
the  real  structure  of  the  attractor  (16),  therefore, 
dimension  estimated  for  an  attractor  of  unnecessarily  high  n 
may  undergo  an  overestimate  since  noise  activity  appearing 
in  the  very  high  dimensional  space  will  get  involved  more  in 
dimensional  estimation. 

Table  4-2  lists  the  correlation  dimension  estimated  for 
four  segments  categorized  as  the  same  brain-state  EEG.  Time 
delay  ranges  from  5 to  11  samples  in  the  case  of  sleep  stage 
0 EEG  segments.  The  variation  is  caused  by  the  inconsistent 
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Table  4-2 

Correlation  Dimension  for  Different  Segments 
Categorized  as  the  Same  Brain  Status 


signal 

segla 

seg2 

seg3 

seg4 

SPIKE (R) 

3.544 

3.408 

3.268 

3.256 

(N=3000) 

(8)b 

(8) 

(9) 

(10) 

SLEEPO(F) 

4.828 

4.835 

4.574 

4.646 

(N=4000) 

(6) 

(5) 

(8) 

(11) 

SLEEPII (R) 

4.770 

4 . 678 

4.705 

4.825 

(N=4000) 

(9) 

(12) 

(12) 

(11) 

SLEEPIII (F) 

3.982 

4.001 

4.034 

3.951 

(N=4000) 

(10) 

(10) 

(12) 

(10) 

SLEEPV(F) 

4.375 

4.294 

3.953 

3.986 

(N=4000) 

(8) 

(8) 

(10) 

(10) 

MUSCLE (R) 

6.952 

6.774 

7.759 

6.651 

(N=4000) 

(6) 

(4) 

(4) 

(7) 

AWSPW(R) 

5.459 

5.560 

5.541 

3 . 119 

(N=4  000 ) 

(17) 

an 

an 

an 

a)  segn:  the  nth  segment. 

b)  number  in  parentheses  shows  the  delay. 
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waveform  patterns  of  the  signal.  Actually,  this  condition 
(the  inconsistency  of  waveform  patterns)  is  also  observed  in 
the  other  brain-state  EEGs  (such  as  the  EEG  segments 
containing  muscle  activity  and  the  awake  EEG  with  spike-and- 
wave) , which  results  in  the  slight  deviation  of  estimated 
correlation  dimension.  This  may  indicate  the  change  of 
brain  dynamics  which  is  not  detected  visually  from  the  one- 
dimensional time  series.  However  the  variance  is  very  small 
for  segments  having  consistent  waveforms.  Note  that  there 
is  a significant  drop  of  the  correlation  dimension  of  the 
fourth  segment  of  awake  EEG  with  spike-and-wave.  It  is 
because  the  spike  has  dominated  the  entire  segment  (Figure 
4-16)  during  that  time  period. 

Largest  Lyapunov  Exponent 

In  our  study  thus  far,  the  estimation  of  the  largest 
Lyapunov  exponent  for  the  EEG  segments  still  undergoes  the 
difficulties  of  determining  the  EVOLV  due  to  the  extremely 
complex  structure  of  the  EEG  attractors  caused  by  a large 
number  of  frequency  components.  Further  investigation  is 
necessary  to  ensure  the  reliability  and  significance  of  the 
estimation.  However,  the  approach  proposed  in  Chapter  3 is 
a potential  way  to  start  with. 

The  procedure  described  in  the  conclusion  of  Chapter  3 
is  repeated  below.  First  the  EVOLV  is  determined  from 
equation  (3.6).  Then  the  spatial  extent  (extent)  is 
computed  as  the  largest  distance  between  any  pair  of  n- 
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dimensional  points  in  the  reconstructed  trajectory.  The 
minimum  scale  limitation  (SCALMN  and  pscalmn)  can  be  derived 
from  the  information  (knee)  of  CIM  curves.  For  each  EEG 
segment,  we  estimate  the  largest  Lyapunov  exponent  for 
pscalmx  ranging  from  10%  to  30%  since  from  our  experiment  we 
find  that  this  pscalmx  range  always  covers  the  range  of 
stationary  estimation.  The  final  result  will  be  computed  by 
averaging  the  exponents  within  the  range  providing 
stationary  estimation.  Table  4-3  lists  the  estimated 
results.  The  range  of  pscalmx  giving  stationary  estimation 
of  Tj  is  shown  in  the  Table.  Note  that  1^  estimated  for 
spike  segment  is  comparably  large.  It  may  be  caused  by  the 
mixed  structure  of  the  background  EEG  activity  and  the 
cyclic  pattern  of  spikes  which  occupy  a very  large  area  in 
phase  space.  The  combination  of  two  different  attractor's 
structure  is  very  likely  to  provide  an  enormous  contribution 
to  the  estimation. 

Discussion 

From  the  observation  in  this  Chapter,  a few  points  may 
be  helpful  to  the  study  of  EEG's  nonlinear  dynamics.  First 
of  all,  we  have  shown  that  the  quantification  of  EEG 
dynamics  through  the  estimation  of  correlation  dimension  has 
potential  to  classify  the  EEG  signals.  From  the  study  of 
estimating  correlation  dimension  for  four  consecutive 
segments,  we  found  that  this  method  can  detect  the  variation 
of  waveform  activity  (such  as  the  case  of  awake  EEG  with 
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Table  4-3 

The  Estimated  Largest  Lyapunov  Exponent 
for  Different  Brain-State  EEG  Segments 


signal 

N 

na 

Tb 

EVOLV 

pscalmn 

pscalmx 

ri 

SPIKE 

3000 

5 

8 

38 

4.6% 

10%-25% 

4 . 20±0 . 34c 

SLEEPO 

4000 

7 

6 

24 

9.3% 

18%-26% 

1 . 84±0 . 48 

SLEEPII 

4000 

6 

9 

51 

7.8% 

14%-24% 

2 . 08±0 . 14 

SLEEPIII 

4000 

5 

10 

54 

4.7% 

12%-20% 

2 . 7 0±0 . 05 

SLEEPV 

4000 

6 

8 

31 

7.2% 

16%-24% 

2 . 50±0 . 34 

MUSCLE 

4000 

11 

6 

9 

4.7% 

24%-30% 

0.46+0.07 

AWSPW 

4000 

8 

17 

60 

6.2% 

14%-2  0% 

0 . 73±0 . 04 

a)  n:  embedding  dimension. 

b)  r:  delay  in  samples. 

c)  The  average  and  deviation  are  computed  from  the  exponents 
estimated  from  the  specified  pscalmx  range. 
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spike-and-wave,  Table  4-2) . In  addition,  we  have 
corroborated  that  the  information  obtained  from  the  signal 
(such  as  the  power  spectrum)  may  be  used  to  understand  the 
characteristics  of  its  reconstructed  trajectory.  The 
quantified  topological  properties  of  the  reconstructed 
trajectory  provides  a helpful  measurement  for  the  change  of 
activity  of  the  original  one-dimensional  time  series.  The 
results  of  our  study  in  this  Chapter  show  that  /i  (MUSCLE)  > 

H (AWSPW)  > n (SLEEP  II)  > /x  (SPIKE)  (estimated  from  the  raw 
segments)  and  /i( SLEEP  0)  > n (SLEEP  V)  > n( SLEEP  III) 
(estimated  from  the  filtered  segments) , while  the  reverse 
order  is  observed  for  the  estimated  largest  Lyapunov 
exponent  (AWSPW:  awake  EEG  with  spike-and-wave) . Note  that 
the  AWSPW  segment  selected  here  (the  first  segment)  involves 
little  spike  activity  (Figure  4-7) . Hence  a large  number  of 
correlation  dimension  is  obtained  («  5.5).  One  point  needed 
to  be  aware  is  that  the  estimation  of  the  largest  Lyapunov 
exponent  might  suffer  a questionable  reliability  when  the 
attractor  is  of  very  high  dimension  (16) . However,  we  have 
not  found  any  numerical  evidence  specifying  the  upper  bound 
of  dimension  allowed  when  employing  this  methodology 
(Appendix)  in  experimental  data. 

Although  the  results  are  encouraging,  they  are 
preliminary  since  the  data  set  is  small.  For  the  future 
study,  the  methods  of  dimension  and  exponent  estimate  may  be 
employed  to  analyzing  a large  subject  population  and 
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measuring  the  effect  of  inter  and  intra  subject  variability 
in  the  estimated  correlation  dimension  and  Lyapunov 
exponent . 


CHAPTER  5 

CHAOTIC  DYNAMICS  OF  TIME-DELAY  NEURAL  NETWORKS 
Artificial  neural  networks  have  been  intensively 
studied  and  applied  primarily  for  the  generation,  retrieval, 
storage,  and  recognition  of  steady  patterns.  In  those 
cases,  the  dynamics  of  the  network  always  converge  to  a 
fixed-point  attractor.  This  dynamic  behavior  cannot  be  used 
to  describe  the  EEGs  with  the  well  known  continuous  and 
spontaneously  changing  activities.  In  order  to  develop  an 
artificial  system  to  simulate  the  dynamic  behavior  of  the 
EEG,  we  have  to  either  modify  the  existing  models  or  explore 
a new  model  capable  of  generating  complex  oscillatory 
activity.  For  decades,  biologists  and  neurologists  have 
tried  to  investigate  the  neuronal  basis  of  rhythmic  behavior 
through  the  analysis  and  modelling  of  oscillatory  neural 
networks  in  invertebrates  (84,85).  From  former  experiments 
a principal  conclusion  emerges,  i.e.  the  cellular  and 
synaptic  properties  of  the  neural  networks  are  most 
important  to  the  successful  functioning  of  oscillators 
(cyclic-pattern  generators)  (84) . This  concept  actually  is 
one  of  the  reasons  for  the  artificial  neural  networks 
resurgence  and  is  believed  to  lead  to  the  sixth-generation 
technology.  However,  most  studies  deal  with  the  problem 
(association  and  memory  retrieval)  of  searching  a set  of 
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weights  such  that  the  network  has  fixed-point  attractors  at 
given  points  (86) . Kleinfeld  (33)  , and  Sompolinsky  and 
Kanter  (34)  introduced  a class  of  networks  (time-delay 
recurrent  neural  networks)  capable  of  generating  a 
memorizing  seguential  patterns.  The  system  dynamics 
traverse  phase  space  in  a limit  cycle,  instead  of  resting  in 
an  eguilibrium  state.  Chapter  1 has  described  some 
characteristics  of  the  time-delay  neural  networks  as  well  as 
the  approach  of  constructing  the  networks,  which  will  be 
called  hereinafter  the  "T, D networks"  since  the  model  is 
characterized  by  two  sets  of  weights,  T^j  which  accounts  for 
the  effect  of  the  current  state,  and  Di j which  makes  the 
former  output  state  (Kdly  steps  earlier)  of  a processing 
element  (PE)  a part  of  the  present  activity  (input  state  of 
the  PE) . 

In  this  study,  our  major  goal  is  to  study  the  nonlinear 
dynamical  properties  of  the  T,D  networks.  Understanding  the 
dynamical  characteristics  of  artificial  neural  networks  may 
lead  to  a new  approach  to  study  and  model  brain  dynamics. 

To  our  belief,  if  an  artificial  neural  network  model 
possesses  an  attractor  with  the  characteristics  similar  to 
the  observed  phase  portrait  of  EEG,  then  the  dynamical 
properties  of  the  brain  state  can  be  studied  by  the  patterns 
of  activation  of  the  artificial  network  model.  This  study 
promotes  the  understanding  and  characterization  of  brain 
status  through  signal  modelling  and  simulation. 
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Experimental  Implementations 

To  speed  up  the  network  simulation  by  eliminating  the 
time-consuming  process  of  numerically  solving  the 
differential  equations,  the  simplified  case  described  below 
will  be  applied  to  our  preliminary  investigation  of  the 
networks.  Let  us  assume  the  charging  time  (time  constant  of 
the  PE)  rs  « 0 and  the  delay  function  impulsional,  i.e. 

D (t , Tp)  = S (t-rD)  . The  discrete  equations  describing  the 
network  dynamics  become 

N N 

Ui(k)  = S TijVWk)  + S D^VWk-ko),  (5.1) 

j=l  j=l 

where  the  evolved  step  k ~ rs  and  the  delay  interval  kD  « 
td/ts  (33)  is  supposed  to  be  the  same  for  all  the  PEs.  As 
described  in  Chapter  1,  the  weights  Tij  and  Dij  are 
determined  from  the  embedded  patterns  (or  memories) . The 
output  state  Vj  (k)  in  our  simulation  algorithm  is  evaluated 
by  the  sigmoid  function  (76,77) 

vi  = g(Uj.)  = (2/tt)  tan“1(7raUi/2)  , (5.2) 

where  the  gain  scaling  parameter  a controls  the  steepness  of 
the  sigmoid  curve.  We  will  work  with  a = 1.5  in  our  models. 
In  the  following  study,  we  concentrate  on  the  T, D network 
with  20  processing  elements  (Nneuron  = 20)  and  5 embedded 
patterns  (Npattern  = 5) . 
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First  of  all  the  dynamical  properties  (correlation 
dimension  and  the  largest  Lyapunov  exponent)  of  the  T, D 
networks  with  fixed  weights  (i.e.  constant  weights)  are 
evaluated.  The  fixed  weight  T, D networks  represent  the 
networks  having  constant  and  throughout  the  time 
evolution.  From  the  work  of  Kleinfeld  (33)  and  Sompolinsky 
et  al.  (34),  the  output  sequences  of  this  class  of  model 
neural  networks  show  cyclic  patterns  (i.e.  a limit  cycle) . 

A limit  cycle,  being  a singularity  of  order  one,  means  that 
the  networks  are  "not  chaotic  enough."  The  CIM  curves  can 
provide  a good  illustration  of  this  fact.  In  our  study,  we 
consider  the  overlap  of  the  current  state  with  the  delayed 
state,  <V(k) • V(k-kD) >,  as  the  output  of  the  entire  network 
because  the  overlap  defined  as 


N 

<V (k) • V (k-kD) > = ( 1/N)  2 Vj (k)  Vj  (k-k„) 


(5.3) 


measures  the  integrated  activities  of  all  the  PEs  in  the 
network.  The  correlation  dimension  and  the  largest  Lyapunov 
exponent  are  estimated  for  the  generated  waveform  <V(k)-V(k- 
kD)  > • 

The  delay  interval  kD  = 10  is  used.  The  parameter 
lsequence  (<  Npattern)  for  computing  the  is  set  at 
lsequence  = Npattern  -1=4.  One  point  needed  to  be 
mentioned  is  that  the  set  of  memories  used  to  form  the  T and 
D matrices  will  not  be  restricted  to  the  orthogonal 
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patterns,  since  it  is  not  our  objective  to  accurately 
retrieve  the  stored  patterns.  The  bit  patterns  of  the 
memories  are  therefore  chosen  at  random.  Based  on  our 
previous  study  of  the  dynamical  properties  of  sleep  stage  II 
EEG,  the  15-second  segment  already  provides  a convergent  and 
reliable  estimate.  The  waveforms  generated  by  the  T,D 
network  are  not  supposed  to  be  more  complicated  or  chaotic 
than  the  EEGs.  Accordingly,  the  segment  length  of  3000 
samples  is  sufficient  for  the  dimension  and  exponent 
estimates. 

Network  Simulation  and 
Estimation  of  the  Dynamical  Properties 

Two  sets  of  embedded  patterns  are  plotted  in  Figure  5- 
1(a)  and  5-1 (b) . The  bottom  trace  (Figure  5-1 (c) ) displays 
the  initial  state  of  both  networks  (network  A and  B) 
constructed  from  memory  set  A (Figure  5-1 (a))  and  set  B 
(Figure  5-1 (b) ) . Note  that  only  two  bits  in  set  A and  B are 
(the  first  bit  in  M2  and  the  second  bit  in  M2) . 

From  eguation  (1.13),  the  little  variation  of  two  memory 
sets  actually  causes  an  18.5%  difference  in  number  of 
weights  between  memory  set  A (Figure  5-2)  and  set  B (Figure 
5“3) . Figure  5-2  and  5-3  display  the  weight  distribution 
maps  of  T-weight  and  D-weight  for  network  A and  B 
respectively.  The  jth  row  depicts  the  weight  of  T^  (D^j),  i 

— 1,  ...,  Nneuron.  Diagonal  indicates  the  weight  T^j  or  D^j 

- 0.  Area  of  each  box  is  proportional  to  the  magnitude  of 
the  corresponding  weight.  A negative  (inhibitory)  weight  is 
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Figure  5-1  Bit  patterns. 

(a)  Memory  set  A;  (b)  Memory  set  B; 
(c)  Initial  state. 
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Figure  5-2  Weight  distribution  maps  for  time-delay 
network  A. 

(a)  Weights  map; 

(b)  Weights  Dij  map. 
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Figure  5-3  Weight  distribution  maps  for  time-delay 
network  B. 

(a)  Weights  map; 

(b)  Weights  Dij  map. 
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Figure  5-3  (continued) . 
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described  by  a shaded  box.  Apparently,  a large  portion  of 
D-weight  is  zero,  whereas  only  Ti:L  = 0 in  T-weight. 

In  Figure  5-4  and  5-5,  we  show  the  800-point  waveforms 
of  the  overlap  of  the  network's  current  state  with  each  of 
the  (five)  memories  for  network  A and  B respectively,  i.e., 
<V(k)-MAi>  and  <V(k)-MBi>,  i = 1,  2,  3,  4,  and  5.  The 
overlaps  display  the  periodic  patterns  with  different 
frequencies  for  the  two  networks.  According  to  (33)  , the 
network  is  in  the  ith  memory  when  the  waveform  of  overlap 
<V(k)-M1>  reaches  the  magnitude  of  1.  In  our  case,  the 
waveforms  can  only  oscillate  in  the  range  (-0.8, 0.8)  because 
the  output  function  (equation  (5.2))  is  not  a linear 
threshold  function  (31).  From  equation  (5.2),  a in  the 
nonlinear  sigmoid  function  determines  the  slope  of  the  1-0 
curve  as  well  as  the  maximum  magnitude  of  the  overlap.  Our 
experiment  shows  that  a larger  a may  result  in  a more 
complex  waveform  pattern  in  one  period.  On  the  other  hand, 
if  a is  too  small  (a  <<  1.0),  the  network  does  not  properly 
oscillate,  instead,  only  tiny-amplitude  ripples  are 
observed.  At  this  point,  we  do  not  attempt  to  go  any 
further  about  the  role  of  a in  the  T, D network  since  it  does 
not  make  any  global  alteration  in  network  dynamics.  The 
overlaps  <V(k) • V(k-kD)  > for  both  networks  are  apparently 
periodic  waveforms  (Figure  5-6)  with  about  23  and  40  time 
steps  respectively  in  a period.  Therefore  the  network 
repeats  the  same  output  over  and  over,  corroborating  a 
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Figure  5-4  Overlap  of  the  current  state  with  each  memory  (set  A) . 
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(a)  Network  A;  (b)  Network 
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limit-cycle  attractor  in  state  space.  The  development  of 
parallelism  of  CIM  curves  (Figure  5-7  and  5-8)  obviously 
occurs  in  the  very  beginning  (n  = 2) . The  estimated 
correlation  dimensions  are  1.0950  and  1.1158  for  the 
overlaps  <V(k) • V(k-kD)  > of  network  A and  B (D  = 1 for  a 
limit-cycle  attractor) . The  fractional  dimensions  derived 
indicate  either  that  the  networks  possess  chaotic  property 
or  it  is  caused  by  experimental  error.  After  analyzing  the 
long  train  of  data  points,  we  find  the  numbers  in  one  period 
do  not  exactly  repeat  those  in  the  previous  period,  instead 
they  follow  the  same  trend  of  up-and-down  motions  and 
graphically  appear  to  be  periodic  waveforms.  In  either 
case,  the  network  exhibiting  this  low  dimensionality  is 
still  too  simple  to  describe  the  EEG  dynamics.  If  a system 
is  known  to  be  chaotic,  how  chaotic  it  is  turns  out  to  be 
the  next  point  interesting  us.  In  our  study,  the 
dimensionality  and  the  largest  Lyapunov  exponent  are  used  to 
measure  the  "degree  of  chaos  or  strangeness"  (31) . Our 
following  question  is  whether  we  can  break  through  the 
limitation  (in  terms  of  correlation  dynamics)  of  a given 
network  without  changing  the  architecture.  Considering  an 
approach  as  simple  and  straightforward  as  possible,  we  will 
be  investigating  the  possibility  of  pushing  the  network's 
dynamics  into  a more  chaotic  state  by  changing  some  of  the 
network  weights. 
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Figure  5-7  CIM  curves  for  the  overlap  sequence  in  Figure 
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Figure  5-8  CIM  curves  for  the  overlap  sequence  in  Figure  5-6 (b). 
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Imagining  that  one  (or  some)  weight  in  a T, D network 
happens  to  change  its  strength  and  even  its  excitatory  or 
inhibitory  characteristic,  what  can  we  expect  for  the 
output?  If  the  variation  of  Ti;j  or/and  only  occurs  for 
a finite  number  of  times,  the  network's  dynamics  eventually 
will  return  to  its  old  path,  i.e.  limit-cycle  attractor. 
However  the  orbit  most  likely  will  not  be  what  has  been  in 
the  past.  The  network  continues  its  plain,  repetitive  path 
without  remembering  the  disturbance  ever  added  to  it. 
Therefore,  in  the  long  term  the  network's  dynamics  still 
turns  out  a low-dimensional  attractor  (n  « 1.0). 

The  modeling  and  simulation  of  any  kind  of  biological 
system  behavior  or  activity  take  numerous  attempts.  While 
struggling  towards  our  final  goal,  we  might  develop  the  idea 
of  linking  the  recorded  brain  wave  dynamics  with  the 
mechanism  generating  the  complex  EEG.  The  way  we  are 
looking  at  the  T,D  network  is  not  from  the  point  of  view  of 
information  processing  capabilities  (33)  (a  network  with  its 
weights  determined  from  a given  memory  set  generates  the 
seguential  patterns  which  correspond  to  the  transitions 
between  selected  memories) . The  process  of  weight 
modification  is  rather  a continuous  adjustment  of  the 
weights,  differing  from  those  fixed  weight  models  applied 
for  the  pattern  association.  We  attempt  to  investigate 
three  cases  of  weight  variation  (three  rules  of  activation) : 
(1)  once  T^  and  (i,  j = 1,  ...,  Nneuron)  are  determined 
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from  a set  of  Npattern  embedded  patterns,  one  weight  (Tij  or 
Dij , where  i = iv  and  j = jv)  alternates  periodically  (with 
a period  < the  cycling  period  of  the  network  output)  between 
two  levels,  Wold  and  Wnew  (rule  one  activation)  , (2)  same  as 

rule  one  activation  except  that  one  weight  varies 
periodically  but  subject  to  random  amplitude  (rule  two 
activation,  and  (3)  two  sets  of  weights,  (T^^D^j)  and 
(TiBj,DiBj),  are  derived  from  two  memory  sets  (memory  set  A 
and  B,  Figure  5-1) , then  the  T , D network  evolves  with  its 
weights  varying  periodically  between  (TiAj,DiAj)  and 
(TiBj,DiBj)  (rule  three  activation).  We  use  the  rule  of 
periodic  variation  for  simplicity.  Much  remains  unknown 
about  the  complicated  information  processing  system  of  human 
brain.  Even  if  synaptic  modification  underlies  learning, 
thinking,  memorizing,  . . . , and  other  mental  process  of  the 
brain,  how  synapses  adjust  themselves  is  far  beyond  our 
comprehension.  However,  the  existence  of  clocks  in  the  deep 
brain  structures  is  well  known,  which  reinforces  our  choice 
of  periodic  variation  of  network  weights.  In  each  case,  we 
evaluate  the  correlation  dimension  and  the  largest  Lyapunov 
exponent  of  the  network's  output  signal.  There  will  be  no 
external  input  to  the  network  so  that  we  can  totally 
concentrate  on  the  contribution  of  T and  D,  and  the 
complications  of  the  network  can  be  reduced. 
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Rule  One  Activation 

Figure  5-9  shows  the  circuit  diagram  (modified  from 
(33))  of  a 20-PE  T, D network.  Each  diamond  includes  a 
switch  and  two  resistances,  (diAj,diBj)  for  D-diamond  group 
and  (tiAj,tiBj)  for  T-diamond  group.  and  Dij  in  equation 

(5.1)  are  related  to  the  resistances  t^  and  dij  by  | Ti j | = 
Ri/tij  and  |Dij|  = Ri/dij , respectively,  where  Ri  is  derived 
from 


N N 

1/Ri.  = l/rL  + S 1/tij  + S 1/d^.  (5.4) 

j=l  j=l 

The  output  is  proportional  to  the  overlap  <V(k) • V(k-kD) > 
(Nneuron  times  of  the  overlap) . First  of  all,  the  network 
is  constructed  with  weights  determined  from  memory  set  A 
(Figure  5-l(a)).  Suppose  that  T1  2 varies  periodically,  T2>1 
surely  will  change  together  with  T1>2  since  the  weight 
matrix  [Ti;j]  is  kept  symmetric  (Tij  = T.^)  . Figure  5-6 
offers  a good  criterion  to  determine  the  variation  period  of 
T12  and  T2>1.  Since  we  want  to  disturb  the  repetition  of 
the  same  waveform  pattern,  the  variation  period  of  the 
weight  (=  20  steps)  selected  should  be  smaller  than  the 
number  of  steps  (or  samples)  in  one  cycle  (23  steps) . We 
shall  begin  with  the  simple  case  of  "square  wave  variation", 
i.e.,  T12  and  T2>1  alternate  between  two  weight  values  (Wold 
and  Wnew)  . Starting  from  a Wnew  « Wold  and  gradually 
increasing  the  difference  between  Wnew  and  Wold,  we 
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Figure  5-9  Circuit  diagram  for  the  time-delay  artificial  neural  network. 
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experiment  with  several  values  of  Wnew  to  make  the  network's 
output  exhibit  different  waveform  patterns.  With  noisy 
weights  (i.e.,  tolerable  increment  or  decrement  in  finite 
number  of  and  D£j) , the  cyclic  behavior  of  the  network 
essentially  remains  unperturbed  (33).  Hence  |Wnew  - Wold| 
must  be  large  enough  to  keep  the  network  from  generating  the 
cyclic  seguence.  Wnew  = Wold  + (-20)Wold  is  tested  to  be 
strong  enough  to  drive  the  network  to  local  chaos  (explained 
below) . 

Figure  5-10  displays  the  3090-sample  output  <V(k)*V(k- 
kD)>  with  T12  and  T2>1  alternating  between  -0.05  (Wold)  and 
0.95  (Wnew)  every  20  steps  (varying  period  = 20  time  steps). 
The  output  resembles  EEG  waveforms.  However  at  a close 
scrutiny,  the  output  signal  appears  periodic  with  period  « 

60  samples,  yet  the  waveform  in  one  period  resembles  the 
EEG.  Note  that  the  <V(k) • V(k-kD)  > sequence,  although 
looking  periodic  visually  (Figure  5-10) , never  repeats 
numerically  the  previous  period.  What  are  the  implications 
on  the  dimensional  estimation?  Mutual  information  is  first 
computed  to  determine  the  delay  used  in  estimating  the 
correlation  dimension  for  the  3000-sample  segment.  The 
first  significant  local  minimum  is  located  at  delay  = 4 
samples.  The  plot  of  CIM  curves  for  embedding  dimension  n = 
2 to  15  (Figure  5-11)  is  worthy  of  note.  There  are  two 
distinguishable  portions  that  signify  dual  characteristics 
of  the  network.  In  the  large-scale  domain  (-4.6  < ln[C(r)] 
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Figure  5-10  Output  with  T12  anc*  T2,l  alternating  between  -0.05  and  0.95  (period=20) 
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< -0.4),  the  curves  reveal  the  chaotic  nature  like  the  CIM 
curves  for  the  EEG  segment,  i.e.,  slope  increases  with  n up 
to  some  value  (n  = 8)  and  then  saturates.  Correlation 
dimension  estimated  for  this  region  is  /n  = 3.4763. 

Apparently  the  CIM  curves  in  the  region  between  ln(r)  = 1.9 
and  5.1  exhibit  the  same  features  as  those  curves  for  the 
constant-weight  T, D network  (Figure  5-7  and  5-8) , where 
parallelism  begins  at  n = 2.  A very  small  n (0.5276)  is 
obtained.  It  is  highly  possible  to  associate  this  region 
with  the  periodic-like  characteristic  of  the  waveform. 

Rule  Two  Activation 

For  the  EEG,  the  periodic  activities  observed  in  the 
above  network  never  happen,  except  for  some  short  periodic 
patterns.  The  mental  activities  naturally  involve  intricate 
state-transitions;  and  accordingly,  the  complexity  of  weight 
modification  is  beyond  our  imagination.  In  the  following 
stage,  we  concentrate  on  the  point  how  to  "surpass"  the 
network's  resistance  to  regularity  (periodicity).  Assuming 
that  the  variation  of  synaptic  strength  is  random  in 
magnitude  (instead  of  the  "square-wave  variation" 
investigated  previously)  and  periodic  in  time,  the  output  is 
supposed  to  oscillate  without  repetition  and  the 
dimensionality  of  the  network  is  expected  to  have 
significant  increase.  These  facts  were  experimentally 
observed.  Nevertheless,  randomness  gives  rise  to  an 
undesired  effect  on  the  output  when  compared  to  the  normal 
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EEG.  As  plotted  in  Figure  5-12  (varying  period  =20  time 
steps) , the  signal  is  severely  contaminated  by  the  high- 
frequency  noise.  In  our  algorithm  of  random  number 
generator,  one  parameter  (rdrange  > 0.0)  biases  the  amount 
of  high-frequency  noise.  The  weights  (Tx  2 and  T2,i)  to  be 
varied  are  assumed  to  have  their  strength  confined  in  the 
range  -rdrange  < T 1/2,  T2fl  < rdrange.  The  value  of  rdrange 
used  in  Figure  5-12  is  40.0  which  seems  too  large.  On  the 
other  hand,  the  T,D  network  seems  to  be  unaware  of  the 
variation  of  weights  when  a small  value  of  rdrange  is  set. 
For  instance,  the  waveform  in  Figure  5-13 (a)  (rdrange  = 1.0) 
basically  looks  the  same  as  the  output  of  constant -weight 
T,D  network  A (Figure  5-6(a)).  The  cyclic  patterns  are 
evident.  The  periodic  disturbance  appears  to  be  absorbed 
and  incapable  of  stirring  up  higher  dimensionality  chaotic 
behavior.  While  rdrange  increases,  more  variety  in  waveform 
patterns  is  observed  (Figure  5-13) , yet  the  annoying  high- 
frequency  components  are  introduced  into  the  output. 

The  appropriate  range  of  rdrange  varies  among  all  the 
weights  (Ti;j  and  Di;j);  in  this  case  (T1>2  and  T2>1  vary  with 
period  =20  time  steps)  rdrange  around  3.0  provides  a 
nonperiodic  output  with  less  amount  of  high-frequency  noise. 
Yet  the  waveform  still  shows  the  tendency  to  follow  the 
foregoing  segment  after  a short  period  of  evolution.  From 
our  viewpoint,  an  ideal  time  series  generated  by  the  T,D 
network  should  come  up  with  one  characteristic — chaotic 
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behavior  with  deterministic  nature,  or  the  nonrepetitive 
signal  with  intermittent  token  patterns.  Just  like  what  was 
said  in  Gleick's  book  "Chaos"  (36):  chaos  — deterministic 
and  patterned  — pulls  the  data  into  visible  shapes.  At 
this  point,  it  is  hard  to  picture  in  every  respect  the  high- 
dimensional phase  trajectory.  We  shall  compare  the 
dimension  of  signals  generated  with  a small  rdrange  (3.0) 
and  big  rdrange  (40.0),  from  which  we  may  learn  whether  the 
variation  range  of  synapses  will  cause  any  difference  in 
degree  of  chaos.  However  the  signal  may  need  to  be  filtered 
when  using  a big  rdrange. 

The  output  signals  are  plotted  in  Figure  5-14  and  5-15 
for  rdrange  = 3.0  and  40.0  (filtered  by  Butterworth  lowpass 
filter,  fc  = 60  Hz).  At  this  point,  we  need  to  take  into 
account  the  convergent  capability  of  network's  dynamics 
(chaotic  attractor) , that  is  to  say,  does  a 3000-sample 
segment  provide  convergent  dimensional  estimate?  The 
outputs  in  Figure  5-14  and  5-15  obviously  are  no  more  the 
simple  kind  as  before  (Figure  5-6  and  5-10) . It  is 
necessary  to  find  the  minimum  amount  of  signal  (data 
samples)  required  for  the  completion  of  network's  dynamics, 
because  it  may  take  longer  duration  to  obtain  a reliable 
estimate  for  more  complex  waveforms.  Our  study  shows  that 
the  CIM  curves  for  estimating  longer  output  sequences 
considerably  coincide  with  the  curves  obtained  from  3000 
time  steps.  Accordingly  the  network's  output  evolving  for 


T(l,2)  and  T(2,l)  vary,  period=20,  rdrange=3.0 
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Figure  5-14  Output  with  random  weight  amplitude  (T1>2  & T2  2 , period=20,  rdrange=3.0) 


and  T(21)  vary  ( period  = 20,  rdrange  = 40.0,  fc  = 60  Hz 
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time  step  (k) 

Figure  5-15  Waveform  in  Figure  5-12  filtered  by  Butterworth  LPF  (fc=60  Hz) . 
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3000  steps  is  able  to  provide  a reliable  estimate  of  the 
dimensionality.  From  the  MI  curves,  delay  of  4 and  5 time 
steps  are  selected.  The  CIM  curves  (Figure  5-16  and  5-17) 
for  both  cases  exhibit  some  characteristics  which  are  indeed 
encouraging.  As  the  old  acquaintance  (CIM  for  EEG) , slope 
increases  with  the  embedding  dimension  and  eventually 
reaches  saturation  that  defines  the  correlation  dimension. 
With  rdrange  = 3.0,  the  output  apparently  shows  higher 
dimensionality  (Figure  5-16,  /i  = 5.536  ± 0.094)  than  the 
other  rdrange  (Figure  5-17,  /x  = 4.171  ± 0.077).  Parallelism 
begins  around  n = 9.  There  is  no  significant  knee  on  the 
curves  for  large  n.  Is  the  lower  dimension  of  Figure  5-15 
caused  by  the  filtering?  We  analyze  the  unfiltered  time 
series  (Figure  5-12)  and  find  ju  = 5.932,  slightly  larger 
than  the  one  for  rdrange  = 3.0. 

Regarding  the  stationary  of  exponent  estimate,  it  has 
been  mentioned  in  (16)  that  the  technique  of  estimating  the 
largest  Lyapunov  exponent  from  the  experimental  data  does 
not  work  adequately  for  an  attractor  of  very  high  dimension. 
Since  the  operation  of  the  T, D network  involves  one  factor 
of  randomness  (random  magnitude  of  weights) , the  number  of 
degrees  of  freedom  is  expected  to  be  very  large,  which  may 
become  an  obstruction  to  the  exponent  estimation.  For 
instance,  the  calculated  correlation  exponent  converges  at 
embedding  dimension  n = 9 for  rdrange  = 3.0  and  40.0  with 
filtering,  and  n = 12  for  rdrange  =40.0  without  filtering. 
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Figure  5-16  CIM  curves  for  output  sequence  with  random  weight  amplitude  (rdrange=3) . 
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Figure  5-17  CIM  curves  for  filtered  output  with  random  weight  amplitude  (rdrange=40) . 
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This  high  number  of  embedding  dimension  results  in 
difficulty  of  obtaining  a stationary  estimate. 

Of  course  the  real  biological  system  may  not  have  its 
synapses  adjusted  according  to  such  a plain,  simple  law  like 
the  periodic  variation  in  time  neither  the  random  variation 
in  synaptic  strength.  In  addition,  the  assumption  that  only 
one  pair  of  synapses  varies  is  an  oversimplification. 
Nevertheless  from  the  above  study,  we  may  draw  an  idea  about 
how  a simple  network  designed  to  contain  memories  can 
operate  to  induce  measured  chaotic  behavior  from  a global 
output.  The  next  question  addressed  is  whether  the 
variation  of  different  weights  will  result  in  different 
efficacy  in  network's  dynamics.  One  problem  arises.  It 
would  be  impossible  to  study  the  change  of  dynamical 
properties  caused  by  each  individual  weight  in  the  network. 
There  are  190  pairs  of  T's  and  380  D's  in  a 20-neuron  T, D 
network!  To  avoid  the  tedious  and  repetitious  task,  we 
intend  to  pick  three  pairs  of  T-weight  and  three  D-weight 
covering  different  scales  of  strength.  For  each  D or  T 
pair,  we  execute  the  following  three  steps:  (1)  search  for  a 
good  rdrange  and  simulate  the  network's  operation,  (2)  find 
the  delay  used  for  phase  space  construction,  and  (3) 
evaluate  the  dimensionality  of  the  output  sequence. 

Table  5-1  lists  the  results  of  estimation.  Different 
T-weight  does  not  make  much  difference  in  the  dimensionality 
of  the  output.  Another  feature  observed  in  the  CIM  curves 
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Table  5-1 

/x  and  rx  Estimated  for  Time-Delay  Neural  Networks 
with  Time-Varying  Synapses 


synapse 

rdramge 

delay 

(sample) 

/I  (output) 

ri 

(bits/sec) 

*1,3 
T3 , 1 

3.0 

5 

6.0241 

— 

*1,4 

T4,l 

3.0 

5 

6.4259 

— 

T5,14 

T14,5 

3.0 

5 

6.0059 

— 

D7 , 18 

5.0 

5 

5.3409 

— 

00 

rH 

Q 

40.0 

4 

3.1659 

2 . lla 

D12 , 10 

5.0 

5 

5.6366 

— 

D16 , 12 

5.0 

5 

5.6501 

— 

a)  Parameters  used  for  the  exponent  estimate  are: 

embedding  dimension  = 6,  delay  = 4,  pscalmx  = 25%, 
pscalmn  = 5%,  EVOLV  = 30,  N=3000. 
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of  the  network's  output  is  the  "refusal  to  parallelize" 
phenomenon  which  may  be  caused  by  randomness  included  at  the 
"synaptic  level."  Although  the  correlation  exponent  (Figure 
5-18)  appears  to  be  bounded  by  an  upper  limit  (e.g.  with  n 
increased  up  to  18,  correlation  exponent  < 5.77  when  D12,io 
varies) , it  is  difficult  to  judge  at  which  n the  parallelism 
begins  since  the  slope  of  CIM  does  not  really  saturates,  and 
unlike  the  EEG  segment,  no  upper  knee  appears  on  the  CIM 
curves.  At  this  point,  we  do  not  consider  that  the 
estimation  of  the  largest  Lyapunov  exponent  is  reliable 
because,  besides  the  high  dimensionality  which  makes  the 
technique  impracticable,  we  even  cannot  precisely  determine 
the  number  of  embedding  dimension  of  the  signal.  We  have 
experimented  with  the  algorithm  of  exponent  estimate  and 
find  it  difficult  to  obtain  a stationary  estimation. 

From  Table  5-1,  the  D-weights  do  not  drive  the  network 
farther  away  from  regularity  than  T-weights;  rather  the 
correlation  dimension  /lx (output)  decreases  by  about  0.5. 
Interestingly,  there  is  no  more  contamination  of  high- 
frequency  components  in  the  output  with  rdrange  increasing 
even  up  to  40.0  (Figure  5-19).  The  above  observations 
demonstrate  that  in  the  T, D network  the  D-weight  plays  a 
feebler  role  than  the  T-weight. 

Rule  Three  Activation 

This  is  an  identical  problem  as  the  network  discussed 
in  the  "rule  one  activation",  T12  and  T2>1  alternate  between 


correlation  exponent 
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Figure  5-18  CEM  curves  for  output  sequences  with  random 
weight  amplitude  (period=20,  rdrange=3  for 
Tifj  varying  and  5 for  Difj  varying). 
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two  levels  (Figure  5-10  and  5-11) . However,  now  the  weight 
modification  involves  simultaneously  a large  portion  of  the 
weights.  We  may  still  expect  a periodic-like  output 
waveform.  Suppose  that  two  sets  of  weights  (T^Aj,DiAj)  and 
(Ti  j'Di.Bj)  are  determined  from  memory  set  A (Figure  5-l(a)) 
and  B (Figure  5-l(b)),  the  20-PE  T, D network  in  Figure  5-9 
operates  with  its  weights  alternating  periodically  between 

(Ti  j'Di  j)  and  (Ti.Bj'Di.Bj)f  i.e.,  18.5%  of  the  number  of 
weights  varies  each  time.  Figure  5-20  plots  the  output 
signal  for  varying  period  =20  time  steps.  The  waveform 
does  show  the  tendency  of  repetition,  yet  not  exactly 
duplicating  the  former  segments.  The  correlation  dimension 
is  2.789  ± 0.049,  estimated  with  delay  of  5 time  steps. 
Parallelism  of  CIM  curves  (Figure  5-21)  starts  early  at  n = 
5.  In  addition,  the  CIM  curves  exhibit  guite  the  same 
features  as  those  observed  in  Figure  5-11,  except  that  the 
smooth  portion  (slope  « 0.0)  shrinks. 

It  is  more  practical  to  have  the  weights  varying  among 
those  meaning  (T,D)  sets,  i.e.,  sets  of  weights  representing 
different  stored  memories  or  patterns.  Though  the  result  in 
the  third  case  is  guite  the  same  as  those  observed 
previously  where  only  one  randomly-selected  weight  varies, 
it  might  expand  the  idea  of  linking  patterns  of  activation 
of  neural  networks  to  the  recorded  brain  activity. 

Supposing  that  the  brain  dynamics  involve  a sequential  state 
transition  of  embedded  patterns  Mx  -►  M2  Mn,  the 


synaptic  strength  (T,D)  varies  periodically  (period— 20) 
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Figure  5-20  Output  with  weights  alternating  between  (TiAj,DiAj)  and  (TiBj,DiBj). 
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output  signal  recorded  may  reveal  the  integrated  behavior  of 
the  biological  system  undergoing  the  synaptic  modification 
(T^D1)  -*■  (T2 , D2)  (Tn,Dn)  . Considering  the  periodic- 

like  waveform  during  seizure,  is  it  the  case  that  the  spike- 
and-wave  EEG  during  seizure  is  due  to  the  constant-weight 
brain  dynamics  (such  as  Figure  5-6(a))?  Or,  considering  the 
normal  EEGs , does  the  emergence  of  a typical  token  pattern 
(such  as  the  alpha  or  beta  activity)  correspond  to  the  brain 
dynamics  constituted  by  a specific  sequence  of  embedded 
patterns? 

Conclusion 

In  this  chapter,  we  have  shown  that  the  quantification 
and  characterization  of  the  dynamics  of  the  T, D network  and 
its  output  provide  an  interesting  approach  to  model  and 
investigate  EEG  dynamics.  This  study  shows  that  dynamical 
properties  measured  from  the  ensemble  of  processing  units  is 
related  to  patterns  of  variation  on  the  individual  network 
weights.  In  terms  of  brain  dynamics,  a similar  relation  can 
be  pursued  between  the  EEG  and  patterns  of  activation  in  the 
network  synapses. 

We  concentrated  on  the  networks  having  periodic  change 
of  weights  (corresponding  to  the  synaptic  strengths  of  the 
neural  networks)  since  it  is  known  that  there  are  several 
timing  structures  in  the  brain.  By  changing  the  value  of 
the  weights  in  a periodic  fashion  the  output  of  the  T, D 
networks  displays  correlation  dimension  in  the  EEG  range. 


243 


This  corroborates  the  highly  synchronous  activity  of  the 
brain,  since  in  spite  of  the  huge  number  of  neurons,  their 
dynamics  can  be  simulated  by  a 20-PE  T,D  network. 


CHAPTER  6 
CONCLUSIONS 

Summary  of  Present  Work 

The  investigation  of  EEG  dynamical  properties  by  means 
of  phase  space  construction,  dimensional  analysis,  exponent 
estimation,  etc.  provides  a new  and  viable  approach  for  the 
study  of  brain  activity  as  well  as  simulation  with 
artificial  neural  networks.  From  the  bounded  and  relative 
small  correlation  dimension  and  the  positive  largest 
Lyapunov  exponent  estimated  for  EEG  segments,  we  have  shown 
that  the  brain  waves  could  have  been  produced  by  a low- 
dimensional chaotic  system.  However,  instead  of  simply 
applying  the  existing  techniques  to  the  signal  analysis,  our 
study  began  from  the  investigation  of  the  fundamental 
problems  of  the  implementation  methodology  to  real-world 
experimental  data,  then  led  to  the  establishment  of  better- 
defined  approaches. 

Correlation  Dimension 

We  have  shown  that,  for  many  brain-state  EEG  segments, 
convergent  estimation  of  correlation  dimension  is  possible 
with  less  amount  of  data  points  (while  compared  with  the 
work  of  the  others  (6-9)).  This  helps  us  save  lots  of 
computation  time.  The  procedure  of  dimensional  estimation 
proposed  in  Chapter  2 is  summarized  in  the  following.  For  a 


244 


245 


given  EEG,  correlation  dimension  needs  to  be  estimated  for 
several  N's  (numbers  of  data  points)  to  ensure  the 
reliability  and  convergence  of  estimation.  Time  delay  must 
be  re-evaluated  for  each  segment  having  different  duration. 
By  computing  the  mutual  information  function  (27,59),  an 
appropriate  time  delay  is  selected  as  the  first  local 
minimum  that  is  close  to  the  flat-out  value  of  the  curve. 

CIM  curves  for  embedding  dimension  n ranging  from  2 to  20 
are  shown  to  be  able  to  involve  the  saturating  range  of  CEM 
curve,  which  gives  the  estimation  of  embedding  and 
correlation  dimension  of  EEG.  "Knee"  phenomenon  in  CIM 
curves  was  interpreted  as  the  effect  of  noise  contamination. 
The  appearance  of  parallelism  and  knees  in  CIM  curves  is  the 
primary  indicator  for  us  to  determine  the  dimension  of  the 
signal . 

Noise  component  in  EEG  is  very  likely  to  mislead  the 
estimation  and  the  determination  of  implementation 
parameters  (e.g.  delay) . We  corroborated  that  filtering  had 
the  effect  of  decreasing  correlation  dimension.  Filtering 
is  necessary  to  eliminate  the  effect  of  the  added  noise,  but 
further  research  is  needed  to  correlate  the  effect  of 
filtering  with  dynamical  properties  of  signals.  However, 
the  cutoff  frequency  must  be  selected  carefully  not  to 
distort  the  signal  component. 

Our  study  of  correlation  dimension  of  seven  different 
brain-state  EEGs  (Chapter  4)  shows  that  the  dimensional 
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estimation  has  the  potential  to  classify  the  physiological 
states  of  EEG  signal.  From  our  preliminary  experiment,  the 
order  of  estimated  correlation  dimension  is:  m (MUSCLE)  > 
/x(AWSPW)  > /x  ( SLEEP  II)  > n ( SPIKE)  (original  segments)  and 
/i(SLEEP  0)  > ju ( SLEEP  V)  > m(SLEEP  III)  (filtered  segments). 
Lyapunov  Exponent 

From  our  review  of  the  literature,  no  concrete 
criterion  has  been  found  to  determine  the  implementation 
parameters  for  exponent  estimation.  In  our  study,  a 
systematic  procedure  was  developed  to  help  the  selection  of 
parameters  and  to  quantify  their  effect. 

The  procedure  begins  with  the  selection  of  EVOLV.  The 
value  of  EVOLV  is  highly  related  to  the  orbital  size  of  a 
trajectory.  We  proposed  to  use  the  information  contained  in 
the  power  spectrum  to  study  the  orbital  motion  and  further 
help  us  determine  the  appropriate  EVOLV.  Regarding  the 
value  of  SCALMN  (the  minimum  length  scale  limitation) , the 
published  work  (8)  used  5%  of  extent  (spatial  extent) 
without  giving  any  reason.  We  proposed  to  use  the  knee 
phenomenon  in  CIM  curves  to  compute  SCALMN.  It  is  based  on 
our  interpretation  of  knee  as  the  noise  effect. 

Accordingly,  the  appearance  of  knee  indicates  that  the  noise 
effect  becomes  significant  at  that  length  scale  (SCALMN) . 

It  sets  the  minimum  bound  of  the  initial  separation  of  any 
two  nearby  trajectories.  Unfortunately,  we  are  still  unable 
to  establish  a solid  approach  to  determine  an  adequate 
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SCALMX  (the  maximum  length  scale  limitation) . However, 
convergent  estimation  can  be  obtained  by  estimating  rx  for  a 
wide  range  of  SCALMX.  In  the  study,  we  also  introduced  the 
DPq  map  for  understanding  the  contracting/expanding 
activities  of  local  structures.  This  may  provide  another 
tool  for  selecting  the  parameters  EVOLV  and  SCALMX. 

Network  Simulation 

We  have  developed  the  link  between  the  dynamics  of  the 
complex  biological  system  and  the  controllable  artificial 
neural  network  model,  through  the  analysis  and 
quantification  of  their  dynamical  characteristics.  To  match 
the  dynamic  behavior  of  the  real  biological  system,  the 
concept  of  time-varying  weights  is  introduced  into  the 
network  simulation.  Accordingly,  our  network  model  is  able 
to  produce  higher  dimensional  chaotic  activities,  with  the 
complex  patterns  well  resembling  those  recorded  from  the 
human  scalp  and  the  correlation  dimension  estimated  in  the 
EEG  range.  The  results  are  quite  encouraging  since  we  do 
not  have  to  seek  for  a more  complicated  model  than  a 20- 
neuron  network  to  study  the  brain  dynamics.  Modifying  the 
weights,  but  using  the  constant  topology  and  a small  number 
of  neurons,  the  modified  time-delay  neural  networks  have 
been  shown  to  exhibit  the  dynamical  properties  of  the  EEG. 

Our  study  began  with  the  quantification  of  EEG 
dynamical  properties;  and  hopefully,  it  will  promote  the 
understanding  of  complex  brain  activity  through  network 
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simulation  (a  top-down  procedure).  Freeman  (35,87)  started 
from  the  model  construction,  then  investigated  the  output 
dynamics  to  examine  whether  it  matched  the  dynamic 
characteristics  of  the  biological  neural  system  simulated 
(the  olfactory  system) . His  approach,  from  neuroanatomy  to 
dynamic  behavior,  needs  sufficient  information  of  the 
biological  system  so  that  each  element  in  the  artificial 
network  model  was  "well-founded"  according  to  the  observable 
activity  of  the  biological  system.  However,  with  very 
little  knowledge  of  a much  complex  biological  mechanism 
generating  chaotic  activity,  we  feel  that  it  could  be  more 
practical  to  start  with  the  investigation  of  EEG  dynamics. 
Therefore,  this  study  initiates  our  research  work  of 
"modelling  the  brain  dynamics  by  nonlinear  artificial 
network  simulation."  Some  steps  of  further  investigation, 
which  can  be  immediately  continued  from  this  point,  are 
briefly  described  below. 

Future  Work 

Correlation  Dimension 

Only  a few  EEG  segments  from  the  same  patient  have  been 
evaluated  in  our  preliminary  work.  The  results  are 
encouraging.  However,  the  methodology  developed  in  this 
study  should  be  applied  to  more  EEG  segments  from  different 
subjects  and  physiological  states  to  ensure  its  reliability 
and  to  test  its  capability  in  categorizing  brain  states. 

The  effect  of  linear  filtering  on  the  dynamical  properties 
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must  also  be  further  investigated. 

Lyapunov  Exponents 

As  mentioned  in  Chapter  3 and  4,  the  selection  of 
implementing  parameters  is  still  subject  to  lack  of 
knowledge  of  the  EEG  attractors  which  are  of  very  high 
dimensionality  and  extremely  complicated  orbital  structures. 
Especially  the  EVOLV,  the  determination  still  remains 
uncertain.  Wolf  et  al.  (16)  suggested  that  a wide  range  of 
EVOLV  should  be  checked  for  convergent  exponent  estimates 
when  analyzing  systems  with  unknown  mechanism  for  chaos. 
However,  our  experiment  showed  (Table  3-5)  that  there  are 
several  ranges  of  EVOLV  that  provide  stable  estimates  for 
EEG!  The  guestion  is  then,  which  should  we  pick?  From  the 
signal's  power  spectrum,  scattering  of  the  major  freguency 
components  always  results  in  a large  deviation  of  the 
estimated  (non-stat ionary  estimates)  for  different 

EVOLV ' s . 

Since  the  structure  of  an  EEG  attractor  is  supposed  to 
vary  quite  frequently,  the  attractor  involves  orbits  of 
different  sizes  as  well  as  regions  of  varying  phase  space 
velocities.  Therefore,  a variable  EVOLV  algorithm  may 
provide  another  possibility  for  solving  the  problem  of 
exponent  stability.  The  algorithm  can  be  developed  based  on 
two  criteria:  (1)  reducing  the  space  orientation  change  (see 
definition  in  Appendix)  and  (2)  tracing  the  DPQ  activities 
to  avoid  the  stretchings  and  foldings  as  well  as  restrict 
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the  DPq  to  growing  over  large  (»  SCALMX)  . In  addition,  the 
plot  of  space  orientation  changes  vs.  evolution  time  (EVOLV) 
may  provide  some  useful  information  for  studying  the 
structure  of  EEG  trajectories. 

Artificial  Network  Models 

For  the  time-delay  neural  network  model,  several 
controllable  factors  can  be  adjusted  to  investigate  the 
capabilities  and  performance  of  a T, D network,  such  as  the 
number  of  neurons,  stored  memories  which  determine  the 
weights,  patterns  of  activation  (i.e.,  how  the  weights  are 
varied  in  time) , etc.  Here  the  performance  indicates  the 
waveform  patterns  and  the  dynamical  characteristics  of 
network ' s output . 

It  is  corroborated  in  our  preliminary  study  that  the 
time-delay  neural  network  model  has  the  potential  to  model 
nonlinear  dynamical  behavior  of  the  EEG  signals.  However, 
there  have  emerged  some  other  artificial  network  models, 
such  as  the  Hopfield's  neural  network  model  with  time- 
varying  threshold  function  Gj^  (72)  , which  have  the 
capability  of  generating  the  output  with  deterministic 
dynamics  and  chaotic  nature.  Each  unique  model  captures 
some  particular  characteristics  of  the  real  biological 
behavior.  They  offer  the  alternative  for  our  future  study. 


APPENDIX 

ALGORITHM  FOR  ESTIMATING 
THE  LARGEST  LYAPUNOV  EXPONENT 

Basically  the  algorithm  presented  in  (16)  is  applied  to 

our  study.  Programs  and  detailed  documentation  are 

available  in  the  Computational  Neuroengineering  Laboratory. 

With  slight  modification,  Figure  A and  B display  the  flow 

charts  for  the  computation  of  running  exponent  (lyapunov) . 

The  process  of  the  main  program  (Figure  A)  carries  out 

repeated  cycles  of  propagating  and  replacing  the  principal 

axis  vector  (P  (t)  , Q (t ' ) ) , where  P(t)  = point [ind]  is  the 

fiducial  point  and  Q(t')  = point [ind2]  is  the  best  selected 

nearby  point.  Two  variables  (ind  and  ind2)  are  indices  that 

point  to  these  two  points.  The  logarithm  of  the  ratio  of 

final  to  initial  separation, 

ln{D (P (t+EVOLV) , Q (t ' +EVOLV) ) /D(P(t) , Q (t ' ) ) },  updates  a 
running  average  rate  of  orbital  divergence  or  convergence. 

In  the  following  we  explain  in  detail  the  variables,  arrays, 
functions,  and  some  terms  used  in  the  program  (written  in 
Microsoft  C)  and  flow  charts. 

initialization: 

(1)  read  EEG  data, 

(2)  read  the  implementation  parameters,  including 
number  of  data  points  (N) , embedding  dimension  (n) , 
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Figure  A Flow  chart  for  main  program  estimating  the 
largest  Lyapunov  exponent. 
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Figure  B Flow  chart  for  subroutine  finding  a replacement 
for  the  principal  axis  vector. 
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Figure  B (continued) . 
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time  delay  (r) , spatial  extent  (extent) , pscalmx, 

pscalmn,  and  EVOLV, 

(3)  initialize  the  variables: 

ind=0,  index  that  traces  the  fiducial 

trajectory,  starting  from  point [0]  and 
evolving  till  point [total] , 
iterats=0,  total  number  of  propagation  steps, 
ptime=0,  propagation  time  in  second, 
ptime  = EVOLV*dt*iterats, 
sum=0,  sum  of  instantaneous  running  exponent: 

ln(df/di)  / [EV0LV*dt*ln(2) ] , 

where 

dt  = sampling  interval  = 1/sampling  rate, 

di  = initial  separation  D(P(t) ,Q(t ' ) ) , 

df  = final  separation  D (P (t+EVOLV) , Q (t ' +EVOLV) ) , 

extent: 

the  maximum  distance  of  any  pair  of  n-D  data  points, 
extent  = max{distance(i, j) } , i, j=0, ..., total-1  (i=/j), 


EEG [ i ] : 


array  for  EEG  signal  (i  = 0,  ...,  N-l) , 
point [i] : 

n-D  points  in  the  reconstructed  trajectory 
(i  = 0,  . . . , total-1) , 

total: 

total  number  of  points  in  the  reconstructed 
trajectory, 

total  = N - (n-l) *t  - EVOLV, 

ind2 : 


index  that  points  to  a selected  nearby  point  Q(t'), 
distance (m,n) : 


distance  between  point[m]  and  point[n]. 
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dotprod(i,n,m) : 

inner  product  of  vector  (point [ i] , point [n] ) and 
vector  (point [ i] , point [m] ) , 
lyapunov: 

running  exponent, 
output : 

output  file  contains  (ptime, lyapunov)  to  be  plotted  as 
the  Lexp-t  curve  (e.g.  Figure  3-8). 

Note  that  the  values  of  some  variables  are  selected 
according  to  (16).  For  instance,  angmx  is  initialized  as 
0.3  and  the  SCALMX  can  be  relaxed  up  to  5 times  of  its 
initial  setting  (Figure  B) . Whenever  the  final  separation 
of  the  evolved  principal  axis  vector  is  greater  than  SCALMX, 
a replacement  for  point  Q(t'+EV0LV)  (i.e.  point [ind2])  is 
needed.  Figure  B shows  the  flow  chart  for  searching  the 
best  replacement  point  (point [ ind2 ] ) which  must  satisfy  the 
following  requirements:  (1)  distance  between  point [ind2]  and 
point [m]  (the  fiducial  point)  is  not  greater  than  mscalmx 
(the  relaxed  SCALMX)  and  not  less  than  SCALMN,  and  (2)  the 
orientation  change  (i.e.  thmin,  angle  between  the  original 
principal  axis  vector  and  the  replacement)  is  the  smallest 
one  among  all  the  vectors  (point[m] , point [i] ) , where  i = 0, 
...,  total  (i  f m and  i ^ n) . 
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